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

Correct incorrect format of indels

parent f58c56c0
No related branches found
No related tags found
No related merge requests found
...@@ -3,7 +3,7 @@ package nl.lumc.sasc.biopet.tools ...@@ -3,7 +3,7 @@ package nl.lumc.sasc.biopet.tools
import java.io.File import java.io.File
import java.util import java.util
import htsjdk.samtools.reference.FastaSequenceFile import htsjdk.samtools.reference.{ReferenceSequenceFileFactory, ReferenceSequenceFileWalker, ReferenceSequenceFile, FastaSequenceFile}
import htsjdk.variant.variantcontext.{ Allele, GenotypeBuilder, VariantContextBuilder } import htsjdk.variant.variantcontext.{ Allele, GenotypeBuilder, VariantContextBuilder }
import htsjdk.variant.variantcontext.writer.{ VariantContextWriterBuilder, AsyncVariantContextWriter } import htsjdk.variant.variantcontext.writer.{ VariantContextWriterBuilder, AsyncVariantContextWriter }
import htsjdk.variant.vcf._ import htsjdk.variant.vcf._
...@@ -73,12 +73,23 @@ object GensToVcf extends ToolCommand { ...@@ -73,12 +73,23 @@ object GensToVcf extends ToolCommand {
//TODO: Add info fields //TODO: Add info fields
val infoIt = cmdArgs.inputInfo.map(Source.fromFile(_).getLines()) val infoIt = cmdArgs.inputInfo.map(Source.fromFile(_).getLines())
lazy val fastaFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(cmdArgs.referenceFasta, true, true)
for (genotypeLine <- genotypeIt) { for (genotypeLine <- genotypeIt) {
val genotypeValues = genotypeLine.split(" ") val genotypeValues = genotypeLine.split(" ")
val ref = Allele.create(genotypeValues(3), true) val (start, end, ref, alt) = {
val alt = Allele.create(genotypeValues(4)) val start = genotypeValues(2).toInt
val start = genotypeValues(2).toInt if (genotypeValues(4) == "-") {
val end = ref.length - 1 + start val seq = fastaFile.getSubsequenceAt(cmdArgs.contig, start - 1, start + genotypeValues(4).length - 1)
(start - 1, (start + genotypeValues(4).length - 1),
Allele.create(new String(seq.getBases), true), Allele.create(new String(Array(seq.getBases.head))))
} else {
val ref = Allele.create(genotypeValues(3), true)
(start, (ref.length - 1 + start), Allele.create(genotypeValues(3), true), Allele.create(genotypeValues(4)))
}
}
//val start = genotypeValues(2).toInt
//val end = ref.length - 1 + start
val genotypes = samples.toList.zipWithIndex.map { val genotypes = samples.toList.zipWithIndex.map {
case (sampleName, index) => case (sampleName, index) =>
val gps = Array( val gps = Array(
......
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