Skip to content
Snippets Groups Projects
Commit ea7bdf42 authored by Peter van 't Hof's avatar Peter van 't Hof
Browse files

Added contig

parent cfcb5527
No related branches found
No related tags found
No related merge requests found
......@@ -129,7 +129,7 @@ class Snptest(val root: Configurable) extends BiopetCommandLineFunction with Ref
conditional(fullParameterEstimates, "-full_parameter_estimates") +
optional("-method", method) +
optional("-pheno", pheno)
conditional(summaryStatsOnly, "-summary_stats_only") +
conditional(summaryStatsOnly, "-summary_stats_only") +
conditional(covAll, "-cov_all") +
conditional(covAllContinuous, "-cov_all_continuous") +
conditional(covAllDiscrete, "-cov_all_discrete") +
......
......@@ -18,6 +18,7 @@ package nl.lumc.sasc.biopet.extensions.tools
import java.io.File
import nl.lumc.sasc.biopet.core.{Reference, ToolCommandFunction}
import nl.lumc.sasc.biopet.utils.Logging
import nl.lumc.sasc.biopet.utils.config.Configurable
import org.broadinstitute.gatk.utils.commandline.{Input, Output}
......@@ -43,7 +44,7 @@ class SnptestToVcf(val root: Configurable) extends ToolCommandFunction with Refe
override def beforeGraph(): Unit = {
super.beforeGraph()
if (reference == null) reference = referenceFasta()
if (contig == null) throw new IllegalStateException
if (contig == null) Logging.addError("SnptestToVcf is missing a contig")
//if (outputVcf.getName.endsWith(".vcf.gz")) outputFiles :+= new File(outputVcf.getAbsolutePath + ".tbi")
}
......
......@@ -3,9 +3,9 @@ package nl.lumc.sasc.biopet.tools
import java.io.File
import java.util
import htsjdk.samtools.reference.{FastaSequenceFile, ReferenceSequenceFileFactory}
import htsjdk.variant.variantcontext.writer.{AsyncVariantContextWriter, Options, VariantContextWriterBuilder}
import htsjdk.variant.variantcontext.{Allele, GenotypeBuilder, VariantContextBuilder}
import htsjdk.samtools.reference.{ FastaSequenceFile, ReferenceSequenceFileFactory }
import htsjdk.variant.variantcontext.writer.{ AsyncVariantContextWriter, Options, VariantContextWriterBuilder }
import htsjdk.variant.variantcontext.{ Allele, GenotypeBuilder, VariantContextBuilder }
import htsjdk.variant.vcf._
import nl.lumc.sasc.biopet.utils.ToolCommand
......@@ -48,7 +48,7 @@ object SnptestToVcf extends ToolCommand {
infoHeader match {
case Some(header) => parseLines(header, infoIt, cmdArgs)
case _ => logger.info("No header and records found in file")
case _ => logger.info("No header and records found in file")
}
logger.info("Done")
......@@ -59,9 +59,9 @@ object SnptestToVcf extends ToolCommand {
val headerMap = headerKeys.zipWithIndex.toMap
require(headerKeys.size == headerMap.size, "Duplicates header keys found")
val metaLines = new util.HashSet[VCFHeaderLine]()
for (key <- headerKeys if key != "rsid" if key != "chromosome" if key != "position"
if key != "alleleA" if key != "alleleB" if key != "alleleA")
metaLines.add(new VCFInfoHeaderLine(s"ST_$key", 1, VCFHeaderLineType.String, ""))
for (
key <- headerKeys if key != "rsid" if key != "chromosome" if key != "position" if key != "alleleA" if key != "alleleB" if key != "alleleA"
) metaLines.add(new VCFInfoHeaderLine(s"ST_$key", 1, VCFHeaderLineType.String, ""))
val reference = new FastaSequenceFile(cmdArgs.referenceFasta, true)
require(reference.getSequenceDictionary.getSequence(cmdArgs.contig) != null,
......@@ -76,8 +76,9 @@ object SnptestToVcf extends ToolCommand {
.build)
writer.writeHeader(vcfHeader)
val infoKeys = for (key <- headerKeys if key != "rsid" if key != "chromosome" if key != "position"
if key != "alleleA" if key != "alleleB" if key != "alleleA") yield key
val infoKeys = for (
key <- headerKeys if key != "rsid" if key != "chromosome" if key != "position" if key != "alleleA" if key != "alleleB" if key != "alleleA"
) yield key
var counter = 0
for (line <- lineIt if !line.startsWith("#")) {
......@@ -94,7 +95,7 @@ object SnptestToVcf extends ToolCommand {
.stop(end)
.noGenotypes()
val infoBuilder = infoKeys.foldLeft(builder) { case (a,b) => a.attribute("ST_" + b, values(headerMap(b))) }
val infoBuilder = infoKeys.foldLeft(builder) { case (a, b) => a.attribute("ST_" + b, values(headerMap(b))) }
writer.add(builder.id(rsid).make())
......
......@@ -4,10 +4,10 @@ import java.io.File
import java.util
import htsjdk.samtools.reference.FastaSequenceFile
import nl.lumc.sasc.biopet.core.{BiopetQScript, PipelineCommand, Reference}
import nl.lumc.sasc.biopet.core.{ BiopetQScript, PipelineCommand, Reference }
import nl.lumc.sasc.biopet.extensions.Snptest
import nl.lumc.sasc.biopet.extensions.gatk.{CatVariants, SelectVariants}
import nl.lumc.sasc.biopet.extensions.tools.{GensToVcf, SnptestToVcf}
import nl.lumc.sasc.biopet.extensions.gatk.{ CatVariants, SelectVariants }
import nl.lumc.sasc.biopet.extensions.tools.{ GensToVcf, SnptestToVcf }
import nl.lumc.sasc.biopet.pipelines.gwastest.impute.ImputeOutput
import nl.lumc.sasc.biopet.utils.Logging
import nl.lumc.sasc.biopet.utils.config.Configurable
......@@ -107,8 +107,6 @@ class GwasTest(val root: Configurable) extends QScript with BiopetQScript with R
sv.isIntermediate = true
add(sv)
//TODO: snptest
val snptest = new Snptest(this)
snptest.inputGenotypes :+= sv.outputFile
snptest.inputSampleFiles :+= phenotypeFile
......@@ -118,6 +116,7 @@ class GwasTest(val root: Configurable) extends QScript with BiopetQScript with R
val snptestToVcf = new SnptestToVcf(this)
snptestToVcf.inputInfo = snptest.outputFile.get
snptestToVcf.outputVcf = new File(regionDir, s"${region.chr}-${region.start + 1}-${region.end}.snptest.vcf.gz")
snptestToVcf.contig = region.chr
add(snptestToVcf)
region -> snptestToVcf.outputVcf
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment