Example 96 with CloseableIterator
use of htsjdk.samtools.util.CloseableIterator in project jvarkit by lindenb.
the class IndexCovToSV method doWork.
@Override
public int doWork(final List<String> args) {
final IndexCovUtils indexCovUtils = new IndexCovUtils(this.treshold);
final CharSplitter tab = CharSplitter.TAB;
try {
final SAMSequenceDictionary dict = SequenceDictionaryUtils.extractRequired(this.refFile);
final List<String> samples;
SortingCollection<Cell> sorter1 = SortingCollection.newInstance(Cell.class, new CellCodec(), (A, B) -> A.sortSamplePosition(B), this.sortingDelegate.getMaxRecordsInRam(), this.sortingDelegate.getTmpPaths());
sorter1.setDestructiveIteration(true);
try (BufferedReader r = super.openBufferedReader(oneFileOrNull(args))) {
String line = r.readLine();
if (line == null) {
LOG.error("Cannot read first line of input");
return -1;
}
String[] tokens = tab.split(line);
if (tokens.length < 4 || !tokens[0].equals("#chrom") || !tokens[1].equals("start") || !tokens[2].equals("end")) {
LOG.error("bad first line " + line);
return -1;
}
samples = Arrays.asList(tokens).subList(3, tokens.length);
while ((line = r.readLine()) != null) {
if (StringUtil.isBlank(line))
continue;
tokens = tab.split(line);
if (tokens.length != 3 + samples.size()) {
throw new JvarkitException.TokenErrors(samples.size() + 3, tokens);
}
final SAMSequenceRecord ssr = dict.getSequence(tokens[0]);
if (ssr == null) {
LOG.error(JvarkitException.ContigNotFoundInDictionary.getMessage(tokens[0], dict));
return -1;
}
for (int i = 3; i < tokens.length; i++) {
final double score = Double.parseDouble(tokens[i]);
if (indexCovUtils.getType(score).isAmbigous())
continue;
final Cell cell = new Cell();
cell.tid = ssr.getSequenceIndex();
cell.start = Integer.parseInt(tokens[1]);
cell.end = Integer.parseInt(tokens[2]);
cell.sample_id = i - 3;
cell.score = score;
sorter1.add(cell);
}
}
// while readline
sorter1.doneAdding();
}
/* merge adjacent blocks */
SortingCollection<Cell> sorter2 = SortingCollection.newInstance(Cell.class, new CellCodec(), (A, B) -> A.sortPositionSample(B), this.sortingDelegate.getMaxRecordsInRam(), this.sortingDelegate.getTmpPaths());
sorter2.setDestructiveIteration(true);
try (CloseableIterator<Cell> iter1 = sorter1.iterator()) {
try (PeekableIterator<Cell> peeker1 = new PeekableIterator<>(iter1)) {
while (peeker1.hasNext()) {
final Cell first = peeker1.next();
while (peeker1.hasNext()) {
final Cell second = peeker1.peek();
if (first.sample_id != second.sample_id)
break;
if (first.tid != second.tid)
break;
if (first.end + (this.merge_distance + 1) < second.start)
break;
if (!indexCovUtils.getType(first.score).equals(indexCovUtils.getType(second.score)))
break;
// extend first with end of second
first.end = second.end;
// consumme
peeker1.next();
}
sorter2.add(first);
}
// while peeker1
}
}
sorter1.cleanup();
sorter1 = null;
sorter2.doneAdding();
final Set<VCFHeaderLine> metaData = new HashSet<>();
VCFStandardHeaderLines.addStandardFormatLines(metaData, true, VCFConstants.GENOTYPE_KEY, VCFConstants.GENOTYPE_QUALITY_KEY);
VCFStandardHeaderLines.addStandardInfoLines(metaData, true, VCFConstants.END_KEY, VCFConstants.ALLELE_FREQUENCY_KEY, VCFConstants.ALLELE_COUNT_KEY, VCFConstants.ALLELE_NUMBER_KEY);
final VCFFilterHeaderLine filterHomDel = new VCFFilterHeaderLine("HOM_DEL", "There is one Homozygous deletion.");
metaData.add(filterHomDel);
final VCFFilterHeaderLine filterHomDup = new VCFFilterHeaderLine("HOM_DUP", "There is one Homozygous duplication.");
metaData.add(filterHomDup);
/**
* raw value in indexcov
*/
final VCFFormatHeaderLine foldHeader = new VCFFormatHeaderLine("F", 1, VCFHeaderLineType.Float, "Relative number of copy: 0.5 deletion 1 normal 2.0 duplication");
metaData.add(foldHeader);
final VCFInfoHeaderLine infoSvType = new VCFInfoHeaderLine(VCFConstants.SVTYPE, 1, VCFHeaderLineType.String, "SV type");
metaData.add(infoSvType);
final VCFInfoHeaderLine infoSvLen = new VCFInfoHeaderLine("SVLEN", 1, VCFHeaderLineType.Integer, "ABS(SV length)");
metaData.add(infoSvLen);
/*
final VCFInfoHeaderLine infoSvLens = new VCFInfoHeaderLine("SVLENS", 1, VCFHeaderLineType.Integer,"Allele length");
metaData.add(infoSvLens);
final VCFInfoHeaderLine infoSvEnds = new VCFInfoHeaderLine("SVENDS", 1, VCFHeaderLineType.Integer,"Allele ends");
metaData.add(infoSvEnds);
*/
final VCFHeader vcfHeader = new VCFHeader(metaData, samples);
JVarkitVersion.getInstance().addMetaData(this, vcfHeader);
vcfHeader.setSequenceDictionary(dict);
try (VariantContextWriter vcw = this.writingDelegate.dictionary(dict).open(outputFile)) {
vcw.writeHeader(vcfHeader);
final Allele REF_ALLELE = Allele.create("N", true);
try (CloseableIterator<Cell> iter2 = sorter2.iterator()) {
try (EqualRangeIterator<Cell> peeker2 = new EqualRangeIterator<>(iter2, (A, B) -> A.sortChromStart(B))) {
while (peeker2.hasNext()) {
final List<Cell> array0 = peeker2.next();
final Set<Integer> distinct_ends = array0.stream().map(C -> C.end).collect(Collectors.toCollection(TreeSet::new));
final List<List<Cell>> array_of_array;
if (this.normalize_one_allele && distinct_ends.size() > 1) {
array_of_array = new ArrayList<>(distinct_ends.size());
for (final int end_variant : distinct_ends) {
array_of_array.add(array0.stream().filter(C -> C.end == end_variant).collect(Collectors.toList()));
}
} else {
array_of_array = Collections.singletonList(array0);
}
for (final List<Cell> array : array_of_array) {
int an = 0;
if (array.stream().map(C -> indexCovUtils.getType(C.score)).noneMatch(C -> C.isVariant()))
continue;
final Cell first = array.get(0);
final int max_end = array.stream().mapToInt(C -> C.end).max().getAsInt();
final List<Genotype> genotypes = new ArrayList<>(array.size());
final VariantContextBuilder vcb = new VariantContextBuilder();
vcb.chr(dict.getSequence(first.tid).getSequenceName());
vcb.start(first.start);
vcb.stop(max_end);
vcb.attribute(VCFConstants.END_KEY, max_end);
final List<AltAllele> altAlleles = new ArrayList<>();
boolean has_hom_del = false;
boolean has_hom_dup = false;
for (int i = 0; i < array.size(); i++) {
final Cell cell = array.get(i);
final String sampleName = samples.get(cell.sample_id);
an += 2;
final GenotypeBuilder gb;
final IndexCovUtils.SvType type = indexCovUtils.getType(cell.score);
final Allele allele;
AltAllele altAllele;
if (type.isDeletion() || type.isDuplication()) {
altAllele = altAlleles.stream().filter(C -> C.end == cell.end && C.type.equals(type)).findFirst().orElse(null);
if (altAllele == null) {
allele = Allele.create("<" + (type.isDeletion() ? "DEL" : "DUP") + "." + (1 + altAlleles.size()) + ">", false);
altAllele = new AltAllele(type, cell.end, allele);
altAlleles.add(altAllele);
} else {
allele = altAllele.allele;
}
} else {
altAllele = null;
allele = null;
}
switch(type) {
case REF:
gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, REF_ALLELE));
break;
case HET_DEL:
gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, allele));
altAllele.ac++;
break;
case HOM_DEL:
gb = new GenotypeBuilder(sampleName, Arrays.asList(allele, allele));
vcb.filter(filterHomDel.getID());
altAllele.ac += 2;
has_hom_del = true;
break;
case HET_DUP:
gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, allele));
altAllele.ac++;
break;
case HOM_DUP:
gb = new GenotypeBuilder(sampleName, Arrays.asList(allele, allele));
vcb.filter(filterHomDup.getID());
altAllele.ac += 2;
has_hom_dup = true;
break;
default:
gb = new GenotypeBuilder(sampleName, Arrays.asList(Allele.NO_CALL, Allele.NO_CALL));
break;
}
gb.attribute(foldHeader.getID(), cell.score);
if (!type.isAmbigous()) {
gb.GQ(type.getGenotypeQuality(cell.score));
}
genotypes.add(gb.make());
}
// end loop over genotypes
final List<Allele> alleles = altAlleles.stream().map(A -> A.allele).collect(Collectors.toCollection(ArrayList::new));
if (has_hom_del) {
vcb.filter(filterHomDel.getID());
} else if (has_hom_dup) {
vcb.filter(filterHomDup.getID());
} else {
vcb.passFilters();
}
alleles.add(0, REF_ALLELE);
vcb.alleles(alleles);
vcb.genotypes(genotypes);
if (altAlleles.stream().allMatch(C -> C.type.isDeletion())) {
vcb.attribute(VCFConstants.SVTYPE, "DEL");
} else if (altAlleles.stream().allMatch(C -> C.type.isDuplication())) {
vcb.attribute(VCFConstants.SVTYPE, "DUP");
} else {
vcb.attribute(VCFConstants.SVTYPE, "MIXED");
}
vcb.attribute(infoSvLen.getID(), altAlleles.stream().mapToInt(C -> 1 + C.end - first.start).max().orElse(0));
vcb.attribute(VCFConstants.ALLELE_COUNT_KEY, altAlleles.stream().map(A -> A.ac).collect(Collectors.toList()));
vcb.attribute(VCFConstants.ALLELE_NUMBER_KEY, an);
if (an > 0) {
final int final_an = an;
vcb.attribute(VCFConstants.ALLELE_FREQUENCY_KEY, altAlleles.stream().map(A -> A.ac / (double) final_an).collect(Collectors.toList()));
}
vcw.add(vcb.make());
}
}
// end loop while iter
}
}
sorter2.cleanup();
sorter2 = null;
}
return 0;
} catch (final Throwable err) {
LOG.error(err);
return -1;
} finally {
}
}
Also used :
WritingVariantsDelegate(com.github.lindenb.jvarkit.variant.variantcontext.writer.WritingVariantsDelegate)
Genotype(htsjdk.variant.variantcontext.Genotype)
VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine)
DataInputStream(java.io.DataInputStream)
CloseableIterator(htsjdk.samtools.util.CloseableIterator)
Allele(htsjdk.variant.variantcontext.Allele)
Arrays(java.util.Arrays)
CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter)
Program(com.github.lindenb.jvarkit.util.jcommander.Program)
Parameter(com.beust.jcommander.Parameter)
VCFStandardHeaderLines(htsjdk.variant.vcf.VCFStandardHeaderLines)
VCFHeader(htsjdk.variant.vcf.VCFHeader)
TreeSet(java.util.TreeSet)
ParametersDelegate(com.beust.jcommander.ParametersDelegate)
ArrayList(java.util.ArrayList)
HashSet(java.util.HashSet)
DistanceParser(com.github.lindenb.jvarkit.util.bio.DistanceParser)
NoSplitter(com.github.lindenb.jvarkit.util.jcommander.NoSplitter)
DataOutputStream(java.io.DataOutputStream)
StringUtil(htsjdk.samtools.util.StringUtil)
AbstractDataCodec(com.github.lindenb.jvarkit.util.picard.AbstractDataCodec)
Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher)
Path(java.nio.file.Path)
PeekableIterator(htsjdk.samtools.util.PeekableIterator)
VCFConstants(htsjdk.variant.vcf.VCFConstants)
SequenceDictionaryUtils(com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils)
SortingCollection(htsjdk.samtools.util.SortingCollection)
VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine)
VCFHeaderLineType(htsjdk.variant.vcf.VCFHeaderLineType)
GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder)
Logger(com.github.lindenb.jvarkit.util.log.Logger)
SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary)
Set(java.util.Set)
IOException(java.io.IOException)
JVarkitVersion(com.github.lindenb.jvarkit.util.JVarkitVersion)
EOFException(java.io.EOFException)
Collectors(java.util.stream.Collectors)
JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException)
List(java.util.List)
FractionConverter(com.github.lindenb.jvarkit.jcommander.converter.FractionConverter)
EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator)
VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter)
VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine)
VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine)
BufferedReader(java.io.BufferedReader)
SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord)
Collections(java.util.Collections)
VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)
VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine)
ArrayList(java.util.ArrayList)
SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord)
SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary)
EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator)
ArrayList(java.util.ArrayList)
List(java.util.List)
VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine)
VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter)
VCFHeader(htsjdk.variant.vcf.VCFHeader)
HashSet(java.util.HashSet)
VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine)
CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter)
Genotype(htsjdk.variant.variantcontext.Genotype)
GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder)
VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine)
Allele(htsjdk.variant.variantcontext.Allele)
VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)
BufferedReader(java.io.BufferedReader)
PeekableIterator(htsjdk.samtools.util.PeekableIterator)
Example 97 with CloseableIterator
use of htsjdk.samtools.util.CloseableIterator in project jvarkit by lindenb.
the class VCFCombineTwoSnvs method doVcfToVcf.
@Override
protected int doVcfToVcf(final String inputName, File saveAs) {
BufferedReader bufferedReader = null;
htsjdk.variant.variantcontext.writer.VariantContextWriter w = null;
SortingCollection<CombinedMutation> mutations = null;
CloseableIterator<Variant> varIter = null;
CloseableIterator<CombinedMutation> mutIter = null;
final Map<String, SamReader> sample2samReader = new HashMap<>();
PrintWriter bedPeReport = null;
try {
bufferedReader = inputName == null ? IOUtils.openStreamForBufferedReader(stdin()) : IOUtils.openURIForBufferedReading(inputName);
final VCFUtils.CodecAndHeader cah = VCFUtils.parseHeader(bufferedReader);
/* get VCF header */
final VCFHeader header = cah.header;
final List<String> sampleList = header.getSampleNamesInOrder();
this.indexedFastaSequenceFile = ReferenceSequenceFileFactory.getReferenceSequenceFile(this.referencePath);
final SAMSequenceDictionary dict = SequenceDictionaryUtils.extractRequired(this.indexedFastaSequenceFile);
this.rnaSequenceFactory.setContigToGenomicSequence(C -> getGenomicSequenceForContig(C));
if (this.bamIn != null) {
final Set<String> sampleSet = new HashSet<>(sampleList);
/**
* unroll and open bam file
*/
for (final Path bamFile : IOUtils.unrollPaths(Collections.singletonList(this.bamIn.toString()))) {
LOG.info("opening BAM :" + this.bamIn);
final SamReader samReader = SamReaderFactory.makeDefault().referenceSequence(this.referencePath).validationStringency(ValidationStringency.LENIENT).open(this.bamIn);
if (!samReader.hasIndex()) {
samReader.close();
throw new IOException("Sam file is NOT indexed: " + bamFile);
}
final SAMFileHeader samHeader = samReader.getFileHeader();
if (samHeader.getSequenceDictionary() == null || !SequenceUtil.areSequenceDictionariesEqual(dict, samReader.getFileHeader().getSequenceDictionary())) {
samReader.close();
throw new JvarkitException.DictionariesAreNotTheSame(dict, samReader.getFileHeader().getSequenceDictionary());
}
/* get sample name */
String sampleName = null;
for (final SAMReadGroupRecord rg : samHeader.getReadGroups()) {
if (rg.getSample() == null)
continue;
if (sampleName != null && !sampleName.equals(rg.getSample())) {
samReader.close();
throw new IOException(bamFile + " Contains two samples " + sampleName + " " + rg.getSample());
}
sampleName = rg.getSample();
}
if (sampleName == null) {
samReader.close();
LOG.warn("no sample in " + bamFile);
continue;
}
if (!sampleSet.contains(sampleName)) {
samReader.close();
LOG.warn("no sample " + sampleName + " in vcf. Ignoring " + bamFile);
continue;
}
sample2samReader.put(sampleName, samReader);
}
}
loadTranscripts();
this.variants = SortingCollection.newInstance(Variant.class, new VariantCodec(), new VariantComparatorTwo(dict), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
this.variants.setDestructiveIteration(true);
ProgressFactory.Watcher<VariantContext> progress1 = ProgressFactory.newInstance().dictionary(header).logger(LOG).build();
String vcfLine = null;
while ((vcfLine = bufferedReader.readLine()) != null) {
final VariantContext ctx = progress1.apply(cah.codec.decode(vcfLine));
/* discard non SNV variant */
if (!ctx.isVariant() || ctx.isIndel()) {
continue;
}
/* find the overlapping genes : extend the interval of the variant to include the stop codon */
final Collection<Transcript> genes = this.knownGenes.getOverlapping(new Interval(ctx.getContig(), Math.max(1, ctx.getStart() - 3), ctx.getEnd() + 3)).stream().flatMap(L -> L.stream()).collect(Collectors.toList());
final List<Allele> alternateAlleles = ctx.getAlternateAlleles();
/* loop over overlapping genes */
for (final Transcript kg : genes) {
/* loop over available alleles */
for (int allele_idx = 0; allele_idx < alternateAlleles.size(); ++allele_idx) {
final Allele alt = alternateAlleles.get(allele_idx);
challenge(ctx, alt, kg, vcfLine);
}
}
}
progress1.close();
this.variants.doneAdding();
bedPeReport = this.bedPePath == null ? new PrintWriter(new NullOuputStream()) : IOUtils.openPathForPrintWriter(this.bedPePath);
mutations = SortingCollection.newInstance(CombinedMutation.class, new MutationCodec(), new MutationComparatorTwo(dict), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
mutations.setDestructiveIteration(true);
final VCFFilterHeaderLine vcfFilterHeaderLine = new VCFFilterHeaderLine("TwoHaplotypes", "(number of reads carrying both mutation) < (reads carrying variant 1 + reads carrying variant 2) ");
varIter = this.variants.iterator();
@SuppressWarnings("resource") EqualRangeIterator<Variant> eqVarIter = new EqualRangeIterator<>(varIter, new VariantComparatorOne(dict));
ProgressFactory.Watcher<Variant> progress2 = ProgressFactory.newInstance().dictionary(header).logger(LOG).build();
while (eqVarIter.hasNext()) {
final List<Variant> buffer = eqVarIter.next();
if (buffer.size() < 2)
continue;
for (int i = 0; i + 1 < buffer.size(); ++i) {
final Variant v1 = buffer.get(i);
for (int j = i + 1; j < buffer.size(); ++j) {
final Variant v2 = buffer.get(j);
if (v1.codonStart() != v2.codonStart())
continue;
if (v1.positionInCodon() == v2.positionInCodon())
continue;
if (!v1.wildCodon.equals(v2.wildCodon)) {
throw new IllegalStateException();
}
// no sample share the two variants
final Set<Integer> sharedSamplesIdx = v1.getSharedSampleIndexes(v2);
if (sharedSamplesIdx.isEmpty() && !sampleList.isEmpty())
continue;
final StringBuilder combinedCodon = new StringBuilder(v1.wildCodon);
combinedCodon.setCharAt(v1.positionInCodon(), v1.mutCodon.charAt(v1.positionInCodon()));
combinedCodon.setCharAt(v2.positionInCodon(), v2.mutCodon.charAt(v2.positionInCodon()));
final String pwild = PeptideSequence.of(v1.wildCodon).toString();
final String p1 = PeptideSequence.of(v1.mutCodon).toString();
final String p2 = PeptideSequence.of(v2.mutCodon).toString();
final String pCombined = PeptideSequence.of(combinedCodon).toString();
final String combinedSO;
final String combinedType;
/* both AA are synonymous, while combined is not */
if (!pCombined.equals(pwild) && p1.equals(pwild) && p2.equals(pwild)) {
combinedType = "combined_is_nonsynonymous";
if (pCombined.equals("*")) {
/* http://www.sequenceontology.org/browser/current_svn/term/SO:0001587 */
combinedSO = "stop_gained";
} else if (pwild.equals("*")) {
/* http://www.sequenceontology.org/browser/current_svn/term/SO:0002012 */
combinedSO = "stop_lost";
} else {
/* http://www.sequenceontology.org/miso/current_svn/term/SO:0001992 */
combinedSO = "nonsynonymous_variant";
}
} else if (!pCombined.equals(p1) && !pCombined.equals(p2) && !pCombined.equals(pwild)) {
combinedType = "combined_is_new";
if (pCombined.equals("*")) {
/* http://www.sequenceontology.org/browser/current_svn/term/SO:0001587 */
combinedSO = "stop_gained";
} else {
/* http://www.sequenceontology.org/miso/current_svn/term/SO:0001992 */
combinedSO = "nonsynonymous_variant";
}
} else {
combinedType = null;
combinedSO = null;
}
/**
* ok, there is something interesting here ,
* create two new Mutations carrying the
* two variants
*/
if (combinedSO != null) {
/**
* grantham score is max found combined vs (p1/p2/wild)
*/
int grantham_score = GranthamScore.score(pCombined.charAt(0), pwild.charAt(0));
grantham_score = Math.max(grantham_score, GranthamScore.score(pCombined.charAt(0), p1.charAt(0)));
grantham_score = Math.max(grantham_score, GranthamScore.score(pCombined.charAt(0), p2.charAt(0)));
/**
* info that will be displayed in the vcf
*/
final Map<String, Object> info1 = v1.getInfo(v2);
final Map<String, Object> info2 = v2.getInfo(v1);
// filter for this combined: default it fails the filter
String filter = vcfFilterHeaderLine.getID();
final Map<String, Object> combinedMap = new LinkedHashMap<>();
combinedMap.put("CombinedCodon", combinedCodon);
combinedMap.put("CombinedAA", pCombined);
combinedMap.put("CombinedSO", combinedSO);
combinedMap.put("CombinedType", combinedType);
combinedMap.put("GranthamScore", grantham_score);
info1.putAll(combinedMap);
info2.putAll(combinedMap);
final Map<String, CoverageInfo> sample2coverageInfo = new HashMap<>(sample2samReader.size());
final int chromStart = Math.min(v1.genomicPosition1, v2.genomicPosition1);
final int chromEnd = Math.max(v1.genomicPosition1, v2.genomicPosition1);
/* get phasing info for each sample*/
for (final String sampleName : sample2samReader.keySet()) {
final SamReader samReader = sample2samReader.get(sampleName);
final CoverageInfo covInfo = new CoverageInfo();
sample2coverageInfo.put(sampleName, covInfo);
SAMRecordIterator iter = null;
try {
iter = samReader.query(v1.contig, chromStart, chromEnd, false);
while (iter.hasNext()) {
final SAMRecord rec = iter.next();
if (rec.getReadUnmappedFlag())
continue;
if (rec.isSecondaryOrSupplementary())
continue;
if (rec.getDuplicateReadFlag())
continue;
if (rec.getReadFailsVendorQualityCheckFlag())
continue;
// get DEPTh for variant 1
if (rec.getAlignmentStart() <= v1.genomicPosition1 && v1.genomicPosition1 <= rec.getAlignmentEnd()) {
covInfo.depth1++;
}
// get DEPTh for variant 2
if (rec.getAlignmentStart() <= v2.genomicPosition1 && v2.genomicPosition1 <= rec.getAlignmentEnd()) {
covInfo.depth2++;
}
if (rec.getAlignmentEnd() < chromEnd)
continue;
if (rec.getAlignmentStart() > chromStart)
continue;
final Cigar cigar = rec.getCigar();
if (cigar == null)
continue;
final byte[] bases = rec.getReadBases();
if (bases == null)
continue;
int refpos1 = rec.getAlignmentStart();
int readpos = 0;
boolean found_variant1_on_this_read = false;
boolean found_variant2_on_this_read = false;
/**
* loop over cigar
*/
for (final CigarElement ce : cigar.getCigarElements()) {
final CigarOperator op = ce.getOperator();
switch(op) {
case P:
continue;
case S:
case I:
readpos += ce.getLength();
break;
case D:
case N:
refpos1 += ce.getLength();
break;
case H:
continue;
case EQ:
case M:
case X:
for (int x = 0; x < ce.getLength(); ++x) {
if (refpos1 == v1.genomicPosition1 && same(bases[readpos], v1.altAllele)) {
found_variant1_on_this_read = true;
} else if (refpos1 == v2.genomicPosition1 && same(bases[readpos], v2.altAllele)) {
found_variant2_on_this_read = true;
}
refpos1++;
readpos++;
}
break;
default:
throw new IllegalStateException(op.name());
}
/* skip remaining bases after last variant */
if (refpos1 > chromEnd)
break;
}
/* sum-up what we found */
if (found_variant1_on_this_read && found_variant2_on_this_read) {
covInfo.count_reads_having_both_variants++;
} else if (!found_variant1_on_this_read && !found_variant2_on_this_read) {
covInfo.count_reads_having_no_variants++;
} else if (found_variant1_on_this_read) {
covInfo.count_reads_having_variant1++;
} else if (found_variant2_on_this_read) {
covInfo.count_reads_having_variant2++;
}
}
/* end of loop over reads */
} finally {
iter.close();
iter = null;
}
info1.put("N_READS_BOTH_VARIANTS_" + sampleName, covInfo.count_reads_having_both_variants);
info2.put("N_READS_BOTH_VARIANTS_" + sampleName, covInfo.count_reads_having_both_variants);
info1.put("N_READS_NO_VARIANTS_" + sampleName, covInfo.count_reads_having_no_variants);
info2.put("N_READS_NO_VARIANTS_" + sampleName, covInfo.count_reads_having_no_variants);
info1.put("N_READS_TOTAL_" + sampleName, covInfo.count_reads_having_both_variants + covInfo.count_reads_having_no_variants + covInfo.count_reads_having_variant1 + covInfo.count_reads_having_variant2);
info2.put("N_READS_TOTAL_" + sampleName, covInfo.count_reads_having_both_variants + covInfo.count_reads_having_no_variants + covInfo.count_reads_having_variant1 + covInfo.count_reads_having_variant2);
// count for variant 1
info1.put("N_READS_ONLY_1_" + sampleName, covInfo.count_reads_having_variant1);
info1.put("N_READS_ONLY_2_" + sampleName, covInfo.count_reads_having_variant2);
info1.put("DEPTH_1_" + sampleName, covInfo.depth1);
// inverse previous count
info2.put("N_READS_ONLY_1_" + sampleName, covInfo.count_reads_having_variant2);
info2.put("N_READS_ONLY_2_" + sampleName, covInfo.count_reads_having_variant1);
info2.put("DEPTH_2_" + sampleName, covInfo.depth2);
/* number of reads with both variant is greater than
* reads carrying only one variant: reset the filter
*/
if (2 * covInfo.count_reads_having_both_variants > (covInfo.count_reads_having_variant1 + covInfo.count_reads_having_variant2)) {
/* reset filter */
filter = VCFConstants.UNFILTERED;
info1.put("FILTER_1_" + sampleName, ".");
info2.put("FILTER_2_" + sampleName, ".");
} else {
info1.put("FILTER_1_" + sampleName, vcfFilterHeaderLine.getID());
info2.put("FILTER_2_" + sampleName, vcfFilterHeaderLine.getID());
}
}
/* end of loop over bams */
final CombinedMutation m1 = new CombinedMutation();
m1.contig = v1.contig;
m1.genomicPosition1 = v1.genomicPosition1;
m1.id = v1.id;
m1.refAllele = v1.refAllele;
m1.altAllele = v1.altAllele;
m1.vcfLine = v1.vcfLine;
m1.info = mapToString(info1);
m1.filter = filter;
m1.grantham_score = grantham_score;
m1.sampleIndexes.addAll(sharedSamplesIdx);
m1.sorting_id = ID_GENERATOR++;
mutations.add(m1);
final CombinedMutation m2 = new CombinedMutation();
m2.contig = v2.contig;
m2.genomicPosition1 = v2.genomicPosition1;
m2.id = v2.id;
m2.refAllele = v2.refAllele;
m2.altAllele = v2.altAllele;
m2.vcfLine = v2.vcfLine;
m2.info = mapToString(info2);
m2.filter = filter;
m2.grantham_score = grantham_score;
m2.sampleIndexes.addAll(sharedSamplesIdx);
m2.sorting_id = ID_GENERATOR++;
mutations.add(m2);
bedPeReport.print(m1.contig);
bedPeReport.print('\t');
bedPeReport.print(m1.genomicPosition1 - 1);
bedPeReport.print('\t');
bedPeReport.print(m1.genomicPosition1);
bedPeReport.print('\t');
bedPeReport.print(m2.contig);
bedPeReport.print('\t');
bedPeReport.print(m2.genomicPosition1 - 1);
bedPeReport.print('\t');
bedPeReport.print(m2.genomicPosition1);
bedPeReport.print('\t');
// name
bedPeReport.print(v1.transcriptId);
bedPeReport.print('\t');
// score
bedPeReport.print(grantham_score == GranthamScore.getDefaultScore() ? 0 : (int) ((grantham_score / 255.0) * 1000.0));
bedPeReport.print('\t');
final Transcript kg = this.knownGenes.getOverlapping(new Interval(v1.getContig(), v1.genomicPosition1 - 1, v1.genomicPosition1 + 1)).stream().flatMap(L -> L.stream()).filter(P -> P.getContig().equals(v1.contig) && P.getId().equals(v1.transcriptId)).findFirst().orElseThrow(IllegalStateException::new);
// strand1
bedPeReport.print(kg.isNegativeStrand() ? "-" : "+");
bedPeReport.print('\t');
// strand2
bedPeReport.print(kg.isNegativeStrand() ? "-" : "+");
bedPeReport.print('\t');
if (sharedSamplesIdx.isEmpty()) {
bedPeReport.print('.');
} else {
bedPeReport.print(sharedSamplesIdx.stream().map(I -> sampleList.get(I.intValue())).collect(Collectors.joining(";")));
}
bedPeReport.print('\t');
bedPeReport.print(combinedSO);
bedPeReport.print('\t');
bedPeReport.print(String.join(":", pwild, p1, p2, pCombined));
bedPeReport.println();
}
}
}
}
progress2.close();
mutations.doneAdding();
eqVarIter.close();
eqVarIter = null;
varIter.close();
varIter = null;
variants.cleanup();
variants = null;
bedPeReport.flush();
bedPeReport.close();
bedPeReport = null;
final VCFHeader header2 = new VCFHeader(header);
header2.addMetaDataLine(new VCFHeaderLine(getProgramName() + "AboutQUAL", "QUAL is filled with Grantham Score http://www.ncbi.nlm.nih.gov/pubmed/4843792"));
final StringBuilder infoDesc = new StringBuilder("Variant affected by two distinct mutation. Format is defined in the INFO column. ");
final VCFInfoHeaderLine CodonVariantHeader = new VCFInfoHeaderLine("CodonVariant", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, infoDesc.toString());
header2.addMetaDataLine(CodonVariantHeader);
final VCFInfoHeaderLine CodonSampleHeader = new VCFInfoHeaderLine("Samples", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, "Samples that could be affected");
header2.addMetaDataLine(CodonSampleHeader);
JVarkitVersion.getInstance().addMetaData(this, header2);
if (!sample2samReader.isEmpty()) {
header2.addMetaDataLine(vcfFilterHeaderLine);
}
w = this.writingVariantsDelegate.dictionary(dict).open(IOUtil.toPath(saveAs));
w.writeHeader(header2);
ProgressFactory.Watcher<CombinedMutation> progress3 = ProgressFactory.newInstance().dictionary(header).logger(LOG).build();
mutIter = mutations.iterator();
EqualRangeIterator<CombinedMutation> eqRangeMutIter = new EqualRangeIterator<>(mutIter, new MutationComparatorOne(dict));
while (eqRangeMutIter.hasNext()) {
final List<CombinedMutation> mBuffer = eqRangeMutIter.next();
if (mBuffer.isEmpty())
break;
progress3.apply(mBuffer.get(0));
// default grantham score used in QUAL
int grantham_score = -1;
// default filter fails
String filter = vcfFilterHeaderLine.getID();
final CombinedMutation first = mBuffer.get(0);
final Set<String> info = new HashSet<>();
final VariantContext ctx = cah.codec.decode(first.vcfLine);
final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
vcb.chr(first.contig);
vcb.start(first.genomicPosition1);
vcb.stop(first.genomicPosition1 + first.refAllele.length() - 1);
if (!first.id.equals(VCFConstants.EMPTY_ID_FIELD))
vcb.id(first.id);
for (final CombinedMutation m : mBuffer) {
info.add(m.info);
grantham_score = Math.max(grantham_score, m.grantham_score);
if (VCFConstants.UNFILTERED.equals(m.filter)) {
// at least one SNP is ok one this line
filter = null;
}
}
if (!sampleList.isEmpty()) {
vcb.attribute(CodonSampleHeader.getID(), new ArrayList<>(mBuffer.stream().flatMap(S -> S.sampleIndexes.stream()).map(IDX -> sampleList.get(IDX)).collect(Collectors.toSet())));
}
vcb.unfiltered();
if (filter != null && !sample2samReader.isEmpty()) {
vcb.filter(filter);
} else {
vcb.passFilters();
}
vcb.attribute(CodonVariantHeader.getID(), new ArrayList<String>(info));
if (grantham_score > 0) {
vcb.log10PError(grantham_score / -10.0);
} else {
vcb.log10PError(VariantContext.NO_LOG10_PERROR);
}
w.add(vcb.make());
}
progress3.close();
eqRangeMutIter.close();
mutIter.close();
mutations.cleanup();
mutations = null;
return RETURN_OK;
} catch (final Throwable err) {
LOG.error(err);
return -1;
} finally {
CloserUtil.close(this.indexedFastaSequenceFile);
CloserUtil.close(mutIter);
CloserUtil.close(varIter);
CloserUtil.close(bedPeReport);
if (this.variants != null)
this.variants.cleanup();
if (mutations != null)
mutations.cleanup();
this.variants = null;
for (SamReader r : sample2samReader.values()) CloserUtil.close(r);
CloserUtil.close(w);
CloserUtil.close(bufferedReader);
}
}
Also used :
WritingVariantsDelegate(com.github.lindenb.jvarkit.variant.variantcontext.writer.WritingVariantsDelegate)
Allele(htsjdk.variant.variantcontext.Allele)
Program(com.github.lindenb.jvarkit.util.jcommander.Program)
IOUtil(htsjdk.samtools.util.IOUtil)
Transcript(com.github.lindenb.jvarkit.util.bio.structure.Transcript)
VCFHeader(htsjdk.variant.vcf.VCFHeader)
CigarElement(htsjdk.samtools.CigarElement)
CigarOperator(htsjdk.samtools.CigarOperator)
GenomicSequence(com.github.lindenb.jvarkit.util.picard.GenomicSequence)
SAMFileHeader(htsjdk.samtools.SAMFileHeader)
ReferenceSequenceFile(htsjdk.samtools.reference.ReferenceSequenceFile)
DataOutputStream(java.io.DataOutputStream)
AbstractDataCodec(com.github.lindenb.jvarkit.util.picard.AbstractDataCodec)
Map(java.util.Map)
Path(java.nio.file.Path)
CloserUtil(htsjdk.samtools.util.CloserUtil)
PrintWriter(java.io.PrintWriter)
SequenceDictionaryUtils(com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils)
GranthamScore(com.github.lindenb.jvarkit.util.bio.GranthamScore)
IntervalTreeMap(htsjdk.samtools.util.IntervalTreeMap)
SAMRecordIterator(htsjdk.samtools.SAMRecordIterator)
Collection(java.util.Collection)
Logger(com.github.lindenb.jvarkit.util.log.Logger)
Set(java.util.Set)
Collectors(java.util.stream.Collectors)
JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException)
SAMRecord(htsjdk.samtools.SAMRecord)
ReferenceSequenceFileFactory(htsjdk.samtools.reference.ReferenceSequenceFileFactory)
List(java.util.List)
SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord)
VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter)
VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine)
VariantContext(htsjdk.variant.variantcontext.VariantContext)
VCFHeaderLineCount(htsjdk.variant.vcf.VCFHeaderLineCount)
SamReaderFactory(htsjdk.samtools.SamReaderFactory)
VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)
Genotype(htsjdk.variant.variantcontext.Genotype)
VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine)
DataInputStream(java.io.DataInputStream)
VCFUtils(com.github.lindenb.jvarkit.util.vcf.VCFUtils)
Cigar(htsjdk.samtools.Cigar)
CloseableIterator(htsjdk.samtools.util.CloseableIterator)
PeptideSequence(com.github.lindenb.jvarkit.util.bio.structure.PeptideSequence)
SequenceUtil(htsjdk.samtools.util.SequenceUtil)
ContigNameConverter(com.github.lindenb.jvarkit.util.bio.fasta.ContigNameConverter)
VCFIterator(htsjdk.variant.vcf.VCFIterator)
Parameter(com.beust.jcommander.Parameter)
NullOuputStream(com.github.lindenb.jvarkit.io.NullOuputStream)
AcidNucleics(com.github.lindenb.jvarkit.util.bio.AcidNucleics)
HashMap(java.util.HashMap)
OptionalInt(java.util.OptionalInt)
ValidationStringency(htsjdk.samtools.ValidationStringency)
TreeSet(java.util.TreeSet)
ParametersDelegate(com.beust.jcommander.ParametersDelegate)
RNASequenceFactory(com.github.lindenb.jvarkit.util.bio.structure.RNASequenceFactory)
ArrayList(java.util.ArrayList)
HashSet(java.util.HashSet)
LinkedHashMap(java.util.LinkedHashMap)
ContigDictComparator(com.github.lindenb.jvarkit.util.samtools.ContigDictComparator)
Interval(htsjdk.samtools.util.Interval)
DelegateCharSequence(com.github.lindenb.jvarkit.lang.DelegateCharSequence)
IOUtils(com.github.lindenb.jvarkit.io.IOUtils)
Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher)
WeakHashMap(java.util.WeakHashMap)
VCFConstants(htsjdk.variant.vcf.VCFConstants)
Locatable(htsjdk.samtools.util.Locatable)
SortingCollection(htsjdk.samtools.util.SortingCollection)
VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine)
VCFHeaderLineType(htsjdk.variant.vcf.VCFHeaderLineType)
RNASequence(com.github.lindenb.jvarkit.util.bio.structure.RNASequence)
SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary)
ProgressFactory(com.github.lindenb.jvarkit.util.log.ProgressFactory)
IOException(java.io.IOException)
JVarkitVersion(com.github.lindenb.jvarkit.util.JVarkitVersion)
SamReader(htsjdk.samtools.SamReader)
File(java.io.File)
GtfReader(com.github.lindenb.jvarkit.util.bio.structure.GtfReader)
EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator)
BufferedReader(java.io.BufferedReader)
Comparator(java.util.Comparator)
Collections(java.util.Collections)
VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine)
SAMRecordIterator(htsjdk.samtools.SAMRecordIterator)
HashMap(java.util.HashMap)
LinkedHashMap(java.util.LinkedHashMap)
WeakHashMap(java.util.WeakHashMap)
ProgressFactory(com.github.lindenb.jvarkit.util.log.ProgressFactory)
EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator)
LinkedHashMap(java.util.LinkedHashMap)
HashSet(java.util.HashSet)
CigarOperator(htsjdk.samtools.CigarOperator)
CigarElement(htsjdk.samtools.CigarElement)
VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter)
SAMRecord(htsjdk.samtools.SAMRecord)
SAMFileHeader(htsjdk.samtools.SAMFileHeader)
Interval(htsjdk.samtools.util.Interval)
VCFUtils(com.github.lindenb.jvarkit.util.vcf.VCFUtils)
SAMReadGroupRecord(htsjdk.samtools.SAMReadGroupRecord)
VariantContext(htsjdk.variant.variantcontext.VariantContext)
SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary)
SamReader(htsjdk.samtools.SamReader)
NullOuputStream(com.github.lindenb.jvarkit.io.NullOuputStream)
VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine)
VCFHeader(htsjdk.variant.vcf.VCFHeader)
PrintWriter(java.io.PrintWriter)
Path(java.nio.file.Path)
Transcript(com.github.lindenb.jvarkit.util.bio.structure.Transcript)
IOException(java.io.IOException)
VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine)
Allele(htsjdk.variant.variantcontext.Allele)
Cigar(htsjdk.samtools.Cigar)
VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)
BufferedReader(java.io.BufferedReader)
Example 98 with CloseableIterator
use of htsjdk.samtools.util.CloseableIterator in project jvarkit by lindenb.
the class TView method paint.
void paint(final PrintStream out) {
final Colorizer colorizer;
switch(this.formatOut) {
case html:
colorizer = new HtmlColorizer(out);
break;
case tty:
colorizer = new AnsiColorizer(out);
break;
case plain:
colorizer = new Colorizer(out);
break;
default:
throw new IllegalStateException();
}
if (interval == null) {
LOG.warn("No interval defined");
return;
}
final GenomicSequence contigSequence;
final Function<Integer, Character> refPosToBase;
if (indexedFastaSequenceFile != null) {
final SAMSequenceDictionary dict = SAMSequenceDictionaryExtractor.extractDictionary(referenceFile);
if (dict.getSequence(this.interval.getContig()) == null) {
LOG.warn("No interval with contig " + interval + " in REF");
return;
}
contigSequence = new GenomicSequence(indexedFastaSequenceFile, interval.getContig());
refPosToBase = POS -> {
if (POS < 0 || POS >= contigSequence.length())
return 'N';
return contigSequence.charAt(POS);
};
} else {
contigSequence = null;
refPosToBase = POS -> 'N';
}
/**
* test if genomic position is in interval
*/
final Predicate<Integer> testInInterval = new Predicate<Integer>() {
@Override
public boolean test(final Integer pos) {
return interval.getStart() <= pos && pos <= interval.getEnd();
}
};
final int pixelWidth = this.interval.length();
final Map<Integer, Integer> genomicpos2insertlen = new TreeMap<>();
final Map<String, List<SAMRecord>> group2record = new TreeMap<>();
for (final SamReader samReader : this.samReaders) {
SAMRecordIterator iter = samReader.query(this.interval.getContig(), this.interval.getStart(), this.interval.getEnd(), false);
while (iter.hasNext()) {
final SAMRecord rec = iter.next();
if (rec.getReadUnmappedFlag())
continue;
if (rec.getCigar() == null)
continue;
if (getRecordFilter().filterOut(rec))
continue;
if (!rec.getContig().equals(interval.getContig()))
continue;
if (right().apply(rec) < this.interval.getStart())
continue;
if (this.interval.getEnd() < left().apply(rec))
continue;
String group = this.groupBy.getPartion(rec);
if (group == null || group.isEmpty()) {
group = "undefined_" + this.groupBy.name();
}
List<SAMRecord> records = group2record.get(group);
if (records == null) {
records = new ArrayList<>();
group2record.put(group, records);
}
records.add(rec);
// loop over cigar, get the longest insert
int refpos = rec.getAlignmentStart();
for (final CigarElement ce : rec.getCigar().getCigarElements()) {
if (!this.showInsertions)
break;
final CigarOperator op = ce.getOperator();
if (op.equals(CigarOperator.I) && testInInterval.test(refpos)) {
final Integer longestInsert = genomicpos2insertlen.get(refpos);
if (longestInsert == null || longestInsert.compareTo(ce.getLength()) < 0) {
genomicpos2insertlen.put(refpos, ce.getLength());
}
}
if (op.consumesReferenceBases()) {
refpos += ce.getLength();
}
if (refpos > interval.getEnd())
break;
}
}
CloserUtil.close(iter);
CloserUtil.close(samReader);
}
/**
* compute where are the insertions
*/
// LOG.debug(genomicpos2insertlen);
final Predicate<Integer> insertIsPresentAtX = SCREENX -> {
int x = 0;
int ref = interval.getStart();
while (x < pixelWidth) {
if (x > SCREENX)
return false;
final Integer insertLen = genomicpos2insertlen.get(ref);
if (insertLen == null) {
++x;
++ref;
} else {
if (x <= SCREENX && SCREENX < x + insertLen)
return true;
// (+1) I DON'T UNDERSTAND WHY, BUT IT WORKS
x += (insertLen + 1);
++ref;
}
}
return false;
};
final Function<Character, AnsiUtils.AnsiColor> base2ansiColor = BASE -> {
switch(Character.toUpperCase(BASE)) {
case 'A':
return AnsiUtils.AnsiColor.BLUE;
case 'T':
return AnsiUtils.AnsiColor.GREEN;
case 'G':
return AnsiUtils.AnsiColor.CYAN;
case 'C':
return AnsiUtils.AnsiColor.YELLOW;
default:
return null;
}
};
/**
* print interval title
*/
out.println(interval.getContig() + ":" + interval.getStart() + "-" + interval.getEnd());
/**
* paint base position
*/
int ref = this.interval.getStart();
int x = 0;
out.print(margin("POS:"));
while (x < pixelWidth) {
if (insertIsPresentAtX.test(x)) {
colorizer.pen(AnsiUtils.AnsiColor.RED).print("^");
++x;
} else if ((ref - this.interval.getStart()) % 10 == 0) {
final String f = String.format("%d", ref);
for (int i = 0; i < f.length() && x < pixelWidth; ++i) {
colorizer.pen(AnsiUtils.AnsiColor.GREEN).print(f.charAt(i));
if (!insertIsPresentAtX.test(x))
++ref;
++x;
}
} else {
out.print(".");
++ref;
++x;
}
}
out.println();
/* paint ref base */
out.print(margin("REF:"));
ref = this.interval.getStart();
x = 0;
while (x < pixelWidth) {
if (insertIsPresentAtX.test(x)) {
colorizer.paper(AnsiUtils.AnsiColor.YELLOW).print("*");
++x;
} else {
char refBase = refPosToBase.apply(ref - 1);
colorizer.pen(base2ansiColor.apply(refBase)).print(refBase);
++ref;
++x;
}
}
out.println();
/* loop over samples **/
for (final String groupName : group2record.keySet()) {
if (this.maxReadRowPerGroup == 0)
continue;
final ConsensusBuilder consensus = new ConsensusBuilder();
int y_group = 0;
final List<List<SAMRecord>> rows = new ArrayList<>();
out.println(margin(""));
switch(this.layoutReads) {
case name:
{
rows.addAll(group2record.get(groupName).stream().sorted((R1, R2) -> R1.getReadName().compareTo(R2.getReadName())).map(R -> Collections.singletonList(R)).collect(Collectors.toList()));
break;
}
default:
{
/* pileup reads */
for (final SAMRecord rec : group2record.get(groupName)) {
int y = 0;
for (y = 0; y < rows.size(); ++y) {
final List<SAMRecord> row = rows.get(y);
final SAMRecord last = row.get(row.size() - 1);
if (right().apply(last) + this.distance_between_reads < left().apply(rec)) {
row.add(rec);
break;
}
}
if (y == rows.size()) {
final List<SAMRecord> row = new ArrayList<>();
row.add(rec);
rows.add(row);
}
}
break;
}
}
// each row is only one read, so we need to print the groupName
if (layoutReads == LayoutReads.name) {
out.print(margin(groupName));
out.println();
}
/* print each row */
for (final List<SAMRecord> row : rows) {
++y_group;
boolean print_this_line = (this.maxReadRowPerGroup < 0 || y_group <= this.maxReadRowPerGroup);
if (print_this_line) {
// each row is only one read, print the read name
if (layoutReads == LayoutReads.name) {
String readName = row.get(0).getReadName();
if (row.get(0).getReadPairedFlag()) {
if (row.get(0).getFirstOfPairFlag()) {
readName += "/1";
}
if (row.get(0).getSecondOfPairFlag()) {
readName += "/2";
}
}
out.print(margin(readName));
} else {
out.print(margin(y_group == 1 ? groupName : ""));
}
}
ref = interval.getStart();
x = 0;
for (final SAMRecord rec : row) {
int readRef = left().apply(rec);
// pad before record
while (x < pixelWidth && ref < readRef && testInInterval.test(ref)) {
if (!insertIsPresentAtX.test(x))
++ref;
++x;
if (print_this_line)
out.print(' ');
consensus.add(' ');
}
int readpos = 0;
/* get read base function */
final Function<Integer, Character> baseAt = new Function<Integer, Character>() {
@Override
public Character apply(final Integer readpos) {
final byte[] readBases = rec.getReadBases();
if (readBases == SAMRecord.NULL_SEQUENCE)
return 'N';
if (readpos < 0 || readpos >= rec.getReadLength())
return '?';
return (char) readBases[readpos];
}
};
for (final CigarElement ce : rec.getCigar()) {
final CigarOperator op = ce.getOperator();
if (op.equals(CigarOperator.PADDING))
continue;
/* IN INSERTION, only print if showInsertions is true */
if (this.showInsertions && op.equals(CigarOperator.I)) {
int cigarIdx = 0;
while (x < pixelWidth && cigarIdx < ce.getLength()) {
if (testInInterval.test(readRef)) {
final char readbase = baseAt.apply(readpos);
if (print_this_line)
colorizer.paper(AnsiUtils.AnsiColor.RED).print(readbase);
consensus.add(readbase);
++x;
}
++cigarIdx;
++readpos;
}
continue;
}
int cigarIdx = 0;
while (x < pixelWidth && cigarIdx < ce.getLength()) {
colorizer.clear();
// pad before base
while (x < pixelWidth && testInInterval.test(readRef) && (insertIsPresentAtX.test(x))) {
++x;
if (print_this_line)
colorizer.paper(AnsiUtils.AnsiColor.YELLOW).print("*");
consensus.add(' ');
continue;
}
switch(op) {
case I:
{
// if visible, processed above
if (showInsertions)
throw new IllegalStateException();
readpos++;
break;
}
case P:
break;
case H:
{
if (showClip) {
if (testInInterval.test(readRef)) {
if (print_this_line)
colorizer.paper(AnsiUtils.AnsiColor.YELLOW).print('N');
// CLIPPED base not part of consensus
consensus.add(' ');
++x;
}
++readRef;
}
break;
}
case S:
{
if (showClip) {
if (testInInterval.test(readRef)) {
final char readBase = baseAt.apply(readpos);
if (print_this_line)
colorizer.paper(AnsiUtils.AnsiColor.YELLOW).print(readBase);
// CLIPPED base not part of consensus
consensus.add(' ');
++x;
}
++readpos;
++readRef;
} else {
readpos++;
}
break;
}
case D:
case N:
{
if (testInInterval.test(readRef)) {
if (print_this_line)
colorizer.paper(AnsiUtils.AnsiColor.RED).print('-');
// deletion not not part of consensus
consensus.add(' ');
++x;
}
++readRef;
break;
}
case EQ:
case M:
case X:
{
if (testInInterval.test(readRef)) {
final char refBase = Character.toUpperCase(refPosToBase.apply(readRef - 1));
char readBase = Character.toUpperCase(baseAt.apply(readpos));
consensus.add(readBase);
colorizer.pen(base2ansiColor.apply(readBase));
if (op.equals(CigarOperator.X) || (refBase != 'N' && readBase != 'N' && readBase != refBase)) {
colorizer.pen(AnsiUtils.AnsiColor.RED);
} else if (hideBases) {
if (rec.getReadNegativeStrandFlag()) {
readBase = ',';
} else {
readBase = '.';
}
}
if (showReadName) {
final String readName = rec.getReadName();
if (readpos < 0 || readpos >= readName.length()) {
readBase = '_';
} else {
readBase = readName.charAt(readpos);
}
}
if (rec.getReadNegativeStrandFlag()) {
readBase = Character.toLowerCase(readBase);
} else {
readBase = Character.toUpperCase(readBase);
}
if (print_this_line)
colorizer.print(readBase);
++x;
}
++readpos;
++readRef;
break;
}
}
++cigarIdx;
}
}
// end of loop cigar
ref = readRef;
}
// out.println( " "+ref+" "+row.get(0).getAlignmentStart()+" "+row.get(0).getCigarString()+" "+row.get(0).getReadString());
while (x < pixelWidth) {
if (print_this_line)
out.print(" ");
++x;
}
if (print_this_line)
out.println();
consensus.eol();
if (out.checkError())
break;
}
if (out.checkError())
break;
if (!this.hideConsensus && consensus.bases.stream().anyMatch(C -> C.getCoverage() > 0)) {
out.print(margin(groupName + " CONSENSUS"));
x = 0;
ref = interval.getStart();
while (x < consensus.bases.size()) {
final char refBase = Character.toUpperCase(refPosToBase.apply(ref - 1));
final char consensusBase = consensus.bases.get(x).getConsensus();
if (Character.isWhitespace(consensusBase)) {
// nothing
} else if (refBase != 'N' && consensusBase != refBase) {
colorizer.pen(AnsiUtils.AnsiColor.RED);
} else {
colorizer.pen(base2ansiColor.apply(consensusBase));
}
if (!insertIsPresentAtX.test(x))
++ref;
colorizer.print(consensusBase);
++x;
}
out.println();
}
if (this.numCoverageRows > 0) {
int minCov = consensus.bases.stream().mapToInt(C -> C.getCoverage()).min().orElse(0);
final int maxCov = consensus.bases.stream().mapToInt(C -> C.getCoverage()).max().orElse(0);
for (int y = 0; maxCov > 0 && y < this.numCoverageRows; ++y) {
if (minCov == maxCov)
minCov--;
double fract = (maxCov - minCov) / ((double) this.numCoverageRows);
int inverse_y = (this.numCoverageRows - 1) - y;
int d0 = (int) ((fract) * inverse_y);
// int d1 = (int)((fract) * (inverse_y+1));
out.print(margin(y == 0 ? groupName + " " + maxCov : (y + 1 == this.numCoverageRows ? String.valueOf(minCov) : "")));
for (x = 0; x < consensus.bases.size(); ++x) {
int depth = consensus.bases.get(x).getCoverage() - minCov;
colorizer.print(depth >= d0 ? BLACK_SQUARE : ' ');
}
out.println();
}
}
}
if (this.tabixKnownGene != null && this.indexedFastaSequenceFile != null) {
final List<KnownGene> genes = this.tabixKnownGene.getItemsInInterval(this.interval);
if (!genes.isEmpty()) {
out.println(this.knownGeneUri);
for (final KnownGene gene : genes) {
final KnownGene.CodingRNA codingRna = gene.getCodingRNA(contigSequence);
final KnownGene.Peptide peptide = codingRna.getPeptide();
out.print(margin(gene.getName()));
x = 0;
int ref0 = this.interval.getStart() - 1;
while (x < pixelWidth) {
if (insertIsPresentAtX.test(x)) {
out.print("*");
++x;
} else {
char pepChar = ' ';
if (ref0 >= gene.getTxStart() && ref0 < gene.getTxEnd()) {
pepChar = (gene.isPositiveStrand() ? '>' : '<');
int pepIdx = peptide.convertGenomicToPeptideCoordinate(ref0);
if (pepIdx != -1) {
final AminoAcid aa = AminoAcids.getAminoAcidFromOneLetterCode(peptide.charAt(pepIdx));
final String aa3 = aa == null ? "***" : aa.getThreeLettersCode();
final int[] offset = peptide.convertToGenomicCoordinates(pepIdx);
if (offset != null && offset.length == 3 && aa3 != null && aa3.length() == 3) {
if (offset[0] == ref0)
pepChar = aa3.charAt(0);
else if (offset[1] == ref0)
pepChar = aa3.charAt(1);
else if (offset[2] == ref0)
pepChar = aa3.charAt(2);
else
pepChar = '?';
} else {
pepChar = '?';
}
}
}
out.print(pepChar);
++ref0;
++x;
}
}
while (x < pixelWidth) {
out.print(" ");
++x;
}
out.println();
}
}
out.println();
}
/**
* variant section
*/
if (!this.vcfReaders.isEmpty() && !out.checkError()) {
final Function<GenotypeType, Character> gTypeToSymbol = new Function<GenotypeType, Character>() {
@Override
public Character apply(final GenotypeType gt) {
switch(gt) {
case NO_CALL:
return '?';
case HOM_REF:
return '0';
case HET:
return '1';
case HOM_VAR:
return '2';
case MIXED:
return 'm';
case UNAVAILABLE:
return 'u';
default:
return '.';
}
}
};
out.println();
for (final VcfSource r : this.vcfReaders) {
if (out.checkError())
break;
final VCFHeader header = r.vcfFileReader.getHeader();
final CloseableIterator<VariantContext> iter = r.vcfFileReader.query(this.interval.getContig(), interval.getStart(), interval.getEnd());
final List<VariantContext> variants = new ArrayList<>();
while (iter.hasNext()) {
variants.add(iter.next());
}
iter.close();
if (variants.isEmpty())
continue;
out.println(r.vcfFile.toString());
if (header.hasGenotypingData()) {
for (final String sample : header.getSampleNamesInOrder()) {
if (!variants.stream().map(V -> V.getGenotype(sample)).filter(G -> !hideNoCall || (hideNoCall && !G.isNoCall())).filter(G -> !hideHomRef || (hideHomRef && !G.isHomRef())).findAny().isPresent()) {
continue;
}
out.print(margin(sample));
ref = this.interval.getStart();
x = 0;
while (x < pixelWidth) {
if (insertIsPresentAtX.test(x)) {
out.print("*");
++x;
} else {
char refBase = ' ';
for (final VariantContext ctx : variants) {
if (ctx.getStart() == ref) {
final Genotype g = ctx.getGenotype(sample);
if (g.isNoCall() && this.hideNoCall)
continue;
if (g.isHomRef() && this.hideHomRef)
continue;
refBase = gTypeToSymbol.apply(g.getType());
break;
}
}
out.print(refBase);
++ref;
++x;
}
}
out.println();
}
} else // no genotype
{
for (final VariantContext ctx : variants) {
out.print(margin(String.valueOf(ctx.getStart()) + ":" + ctx.getReference().getDisplayString() + "/" + ctx.getAlternateAlleles().stream().map(A -> A.getDisplayString()).collect(Collectors.joining(","))));
ref = this.interval.getStart();
x = 0;
while (x < pixelWidth) {
if (insertIsPresentAtX.test(x)) {
out.print("*");
++x;
} else {
out.print(ctx.getStart() == ref ? '+' : ' ');
++ref;
++x;
}
}
out.println();
}
}
}
}
}
Also used :
VCFHeader(htsjdk.variant.vcf.VCFHeader)
CigarElement(htsjdk.samtools.CigarElement)
CigarOperator(htsjdk.samtools.CigarOperator)
SAMRecordPartition(com.github.lindenb.jvarkit.util.samtools.SAMRecordPartition)
GenomicSequence(com.github.lindenb.jvarkit.util.picard.GenomicSequence)
ReferenceSequenceFile(htsjdk.samtools.reference.ReferenceSequenceFile)
Map(java.util.Map)
XMLStreamException(javax.xml.stream.XMLStreamException)
XMLStreamWriter(javax.xml.stream.XMLStreamWriter)
Path(java.nio.file.Path)
CloserUtil(htsjdk.samtools.util.CloserUtil)
AnsiUtils(com.github.lindenb.jvarkit.ansi.AnsiUtils)
Predicate(java.util.function.Predicate)
SAMRecordIterator(htsjdk.samtools.SAMRecordIterator)
Logger(com.github.lindenb.jvarkit.util.log.Logger)
GenotypeType(htsjdk.variant.variantcontext.GenotypeType)
Set(java.util.Set)
Collectors(java.util.stream.Collectors)
SAMRecord(htsjdk.samtools.SAMRecord)
ReferenceSequenceFileFactory(htsjdk.samtools.reference.ReferenceSequenceFileFactory)
List(java.util.List)
VariantContext(htsjdk.variant.variantcontext.VariantContext)
SamReaderFactory(htsjdk.samtools.SamReaderFactory)
Genotype(htsjdk.variant.variantcontext.Genotype)
AminoAcid(com.github.lindenb.jvarkit.util.bio.AminoAcids.AminoAcid)
CloseableIterator(htsjdk.samtools.util.CloseableIterator)
Parameter(com.beust.jcommander.Parameter)
Function(java.util.function.Function)
ValidationStringency(htsjdk.samtools.ValidationStringency)
ArrayList(java.util.ArrayList)
HashSet(java.util.HashSet)
Interval(htsjdk.samtools.util.Interval)
NoSplitter(com.github.lindenb.jvarkit.util.jcommander.NoSplitter)
AminoAcids(com.github.lindenb.jvarkit.util.bio.AminoAcids)
IOUtils(com.github.lindenb.jvarkit.io.IOUtils)
Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher)
VCFReaderFactory(com.github.lindenb.jvarkit.variant.vcf.VCFReaderFactory)
SAMSequenceDictionaryExtractor(htsjdk.variant.utils.SAMSequenceDictionaryExtractor)
PrintStream(java.io.PrintStream)
Counter(com.github.lindenb.jvarkit.util.Counter)
SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary)
VCFReader(htsjdk.variant.vcf.VCFReader)
IOException(java.io.IOException)
SamReader(htsjdk.samtools.SamReader)
KnownGene(com.github.lindenb.jvarkit.util.ucsc.KnownGene)
TabixKnownGeneFileReader(com.github.lindenb.jvarkit.util.ucsc.TabixKnownGeneFileReader)
XMLOutputFactory(javax.xml.stream.XMLOutputFactory)
SamRecordFilter(htsjdk.samtools.filter.SamRecordFilter)
SamInputResource(htsjdk.samtools.SamInputResource)
TreeMap(java.util.TreeMap)
Closeable(java.io.Closeable)
SamRecordFilterFactory(com.github.lindenb.jvarkit.util.bio.samfilter.SamRecordFilterFactory)
Collections(java.util.Collections)
SAMRecordIterator(htsjdk.samtools.SAMRecordIterator)
ArrayList(java.util.ArrayList)
Predicate(java.util.function.Predicate)
List(java.util.List)
ArrayList(java.util.ArrayList)
Genotype(htsjdk.variant.variantcontext.Genotype)
CigarOperator(htsjdk.samtools.CigarOperator)
CigarElement(htsjdk.samtools.CigarElement)
SAMRecord(htsjdk.samtools.SAMRecord)
GenotypeType(htsjdk.variant.variantcontext.GenotypeType)
KnownGene(com.github.lindenb.jvarkit.util.ucsc.KnownGene)
VariantContext(htsjdk.variant.variantcontext.VariantContext)
SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary)
SamReader(htsjdk.samtools.SamReader)
Function(java.util.function.Function)
VCFHeader(htsjdk.variant.vcf.VCFHeader)
AminoAcid(com.github.lindenb.jvarkit.util.bio.AminoAcids.AminoAcid)
GenomicSequence(com.github.lindenb.jvarkit.util.picard.GenomicSequence)
TreeMap(java.util.TreeMap)
Example 99 with CloseableIterator
use of htsjdk.samtools.util.CloseableIterator in project jvarkit by lindenb.
the class VCFComposite method doVcfToVcf.
protected int doVcfToVcf(final String inputName, final VCFIterator iterin, final VariantContextWriter out) {
final VCFHeader header = iterin.getHeader();
if (!header.hasGenotypingData()) {
LOG.error("No genotypes in " + inputName);
return -1;
}
final GeneExtractorFactory geneExtractorFactory = new GeneExtractorFactory(header);
final List<GeneExtractorFactory.GeneExtractor> extractors = geneExtractorFactory.parse(this.extractorsNames);
if (extractors.isEmpty()) {
LOG.error("no gene extractor found/defined.");
return -1;
}
final Pedigree pedigree;
try {
final Set<String> sampleNames = new HashSet<>(header.getSampleNamesInOrder());
final PedigreeParser pedParser = new PedigreeParser();
pedigree = pedParser.parse(this.pedigreeFile);
if (pedigree == null || pedigree.isEmpty()) {
LOG.error("pedigree missing/empty");
return -1;
}
this.affectedSamples.addAll(pedigree.getAffectedSamples());
this.affectedSamples.removeIf(S -> !sampleNames.contains(S.getId()));
if (this.affectedSamples.isEmpty()) {
LOG.error("No Affected sample in pedigree. " + this.pedigreeFile + "/" + inputName);
return -1;
}
this.unaffectedSamples.addAll(pedigree.getUnaffectedSamples());
this.unaffectedSamples.removeIf(S -> !sampleNames.contains(S.getId()));
if (pedigree.getUnaffectedSamples().isEmpty()) {
LOG.error("No Unaffected sample in " + this.pedigreeFile + "/" + inputName);
return -1;
}
} catch (final IOException err) {
throw new RuntimeIOException(err);
}
header.addMetaDataLine(new VCFInfoHeaderLine(INFO_TAG, VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, "Variant of VCFComposite"));
if (!StringUtils.isBlank(this.filterNonCompositeTag)) {
header.addMetaDataLine(new VCFFilterHeaderLine(this.filterNonCompositeTag, "Not a Variant fir VCFComposite"));
}
final SAMSequenceDictionary dict = header.getSequenceDictionary();
final Comparator<String> contigCmp;
if (dict == null || dict.isEmpty()) {
contigCmp = (A, B) -> A.compareTo(B);
} else {
contigCmp = new ContigDictComparator(dict);
}
final Comparator<VariantContext> ctxComparator = (V1, V2) -> {
int i = contigCmp.compare(V1.getContig(), V2.getContig());
if (i != 0)
return i;
i = Integer.compare(V1.getStart(), V2.getStart());
if (i != 0)
return i;
return V1.getReference().compareTo(V2.getReference());
};
final Comparator<VariantLine> variantLineComparator = (V1, V2) -> {
final int i = ctxComparator.compare(V1.ctx, V2.ctx);
if (i != 0)
return i;
return Long.compare(V1.id, V2.id);
};
long ID_GENERATOR = 0L;
this.vcfDecoder = VCFUtils.createDefaultVCFCodec();
this.vcfDecoder.setVCFHeader(header, VCFHeaderVersion.VCF4_2);
this.vcfEncoder = new VCFEncoder(header, false, true);
SortingCollection<GeneAndVariant> sorting = null;
SortingCollection<VariantLine> outputSorter = null;
try {
LOG.info("reading variants and genes");
/* Gene and variant sorter */
sorting = SortingCollection.newInstance(GeneAndVariant.class, new GeneAndVariantCodec(), GeneAndVariant::compareGeneThenIndex, this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
sorting.setDestructiveIteration(true);
/* Variant sorter */
outputSorter = SortingCollection.newInstance(VariantLine.class, new VariantLineCodec(), variantLineComparator, this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
outputSorter.setDestructiveIteration(true);
/* read input */
while (iterin.hasNext()) {
final VariantContext ctx = iterin.next();
final VariantLine variantLine = new VariantLine(++ID_GENERATOR, ctx);
if (!this.variantJexl.test(ctx)) {
outputSorter.add(variantLine);
continue;
}
if (!acceptVariant(ctx)) {
outputSorter.add(variantLine);
continue;
}
final Set<GeneIdentifier> geneKeys = new HashSet<>();
extractors.stream().map(EX -> EX.apply(ctx)).flatMap(H -> H.keySet().stream()).forEach(KG -> {
geneKeys.add(new GeneIdentifier(KG.getKey(), KG.getGene(), KG.getMethod().replace('/', '_')));
});
if (geneKeys.isEmpty()) {
outputSorter.add(variantLine);
continue;
}
for (final GeneIdentifier gk : geneKeys) {
final GeneAndVariant gav = new GeneAndVariant(gk, variantLine);
gav.gene.contig = ctx.getContig();
sorting.add(gav);
}
}
sorting.doneAdding();
LOG.info("compile per gene");
this.reportGeneWriter = (this.geneReportPath == null ? new PrintWriter(new NullOuputStream()) : IOUtils.openPathForPrintWriter(this.geneReportPath));
this.reportGeneWriter.print("#CHROM");
this.reportGeneWriter.print('\t');
this.reportGeneWriter.print("bed.start");
this.reportGeneWriter.print('\t');
this.reportGeneWriter.print("bed.end");
this.reportGeneWriter.print('\t');
this.reportGeneWriter.print("gene.key");
this.reportGeneWriter.print('\t');
this.reportGeneWriter.print("gene.label");
this.reportGeneWriter.print('\t');
this.reportGeneWriter.print("gene.source");
this.reportGeneWriter.print('\t');
this.reportGeneWriter.print("count.variants");
this.reportGeneWriter.print('\t');
this.reportGeneWriter.print("affected.counts");
this.reportGeneWriter.print('\t');
this.reportGeneWriter.print("affected.total");
this.reportGeneWriter.print('\t');
this.reportGeneWriter.print("affected.samples");
this.reportGeneWriter.println();
this.reportWriter = (this.reportPath == null ? new PrintWriter(new NullOuputStream()) : IOUtils.openPathForPrintWriter(this.reportPath));
this.reportWriter.print("#CHROM");
this.reportWriter.print('\t');
this.reportWriter.print("bed.start");
this.reportWriter.print('\t');
this.reportWriter.print("bed.end");
this.reportWriter.print('\t');
this.reportWriter.print("gene.index");
this.reportWriter.print('\t');
this.reportWriter.print("gene.key");
this.reportWriter.print('\t');
this.reportWriter.print("gene.label");
this.reportWriter.print('\t');
this.reportWriter.print("gene.source");
for (int side = 0; side < 2; ++side) {
this.reportWriter.print('\t');
final String prefix = "variant" + (side + 1) + ".";
this.reportWriter.print(prefix + "start");
this.reportWriter.print('\t');
this.reportWriter.print(prefix + "end");
this.reportWriter.print('\t');
this.reportWriter.print(prefix + "ref");
this.reportWriter.print('\t');
this.reportWriter.print(prefix + "alt");
this.reportWriter.print('\t');
this.reportWriter.print(prefix + "info");
for (final Sample sn : this.affectedSamples) {
this.reportWriter.print('\t');
this.reportWriter.print(prefix + "gt[" + sn.getId() + "].affected");
}
for (final Sample sn : this.unaffectedSamples) {
this.reportWriter.print('\t');
this.reportWriter.print(prefix + "gt[" + sn.getId() + "].unaffected");
}
}
this.reportWriter.println();
// compile data
CloseableIterator<GeneAndVariant> iter2 = sorting.iterator();
EqualRangeIterator<GeneAndVariant> eqiter = new EqualRangeIterator<>(iter2, (A, B) -> A.gene.compareTo(B.gene));
while (eqiter.hasNext()) {
final List<GeneAndVariant> variants = eqiter.next();
scan(variants.get(0).gene, variants.stream().map(L -> L.variant).collect(Collectors.toList()));
for (final GeneAndVariant ga : variants) outputSorter.add(ga.variant);
}
eqiter.close();
iter2.close();
sorting.cleanup();
//
this.reportWriter.flush();
this.reportWriter.close();
this.reportGeneWriter.flush();
this.reportGeneWriter.close();
LOG.info("write variants");
CloseableIterator<VariantLine> iter1 = outputSorter.iterator();
EqualRangeIterator<VariantLine> eqiter1 = new EqualRangeIterator<>(iter1, variantLineComparator);
out.writeHeader(header);
while (eqiter1.hasNext()) {
final List<VariantLine> array = eqiter1.next();
final VariantContext firstCtx = array.get(0).ctx;
final Set<String> set = getAnnotationsForVariant(firstCtx);
final VariantContext outCtx;
final VariantContextBuilder vcb = new VariantContextBuilder(firstCtx);
for (int y = 1; y < array.size(); ++y) {
set.addAll(getAnnotationsForVariant(array.get(y).ctx));
}
if (set.isEmpty()) {
if (StringUtils.isBlank(this.filterNonCompositeTag)) {
// ignore
continue;
} else {
vcb.filter(this.filterNonCompositeTag);
}
} else {
if (!firstCtx.isFiltered()) {
vcb.passFilters();
}
vcb.attribute(INFO_TAG, new ArrayList<>(set));
}
outCtx = vcb.make();
out.add(outCtx);
}
outputSorter.cleanup();
eqiter1.close();
iter1.close();
return 0;
} catch (final Exception err) {
LOG.error(err);
return -1;
}
}
Also used :
WritingVariantsDelegate(com.github.lindenb.jvarkit.variant.variantcontext.writer.WritingVariantsDelegate)
JexlVariantPredicate(com.github.lindenb.jvarkit.util.vcf.JexlVariantPredicate)
Program(com.github.lindenb.jvarkit.util.jcommander.Program)
GeneExtractorFactory(com.github.lindenb.jvarkit.util.vcf.predictions.GeneExtractorFactory)
VCFHeader(htsjdk.variant.vcf.VCFHeader)
VCFEncoder(htsjdk.variant.vcf.VCFEncoder)
VCFHeaderVersion(htsjdk.variant.vcf.VCFHeaderVersion)
DataOutputStream(java.io.DataOutputStream)
StringUtil(htsjdk.samtools.util.StringUtil)
AbstractDataCodec(com.github.lindenb.jvarkit.util.picard.AbstractDataCodec)
Path(java.nio.file.Path)
PrintWriter(java.io.PrintWriter)
AbstractVCFCodec(htsjdk.variant.vcf.AbstractVCFCodec)
Predicate(java.util.function.Predicate)
Logger(com.github.lindenb.jvarkit.util.log.Logger)
Set(java.util.Set)
EOFException(java.io.EOFException)
Collectors(java.util.stream.Collectors)
List(java.util.List)
StringUtils(com.github.lindenb.jvarkit.lang.StringUtils)
VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter)
VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine)
VariantContext(htsjdk.variant.variantcontext.VariantContext)
VCFHeaderLineCount(htsjdk.variant.vcf.VCFHeaderLineCount)
VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)
Genotype(htsjdk.variant.variantcontext.Genotype)
DataInputStream(java.io.DataInputStream)
CloseableIterator(htsjdk.samtools.util.CloseableIterator)
VCFUtils(com.github.lindenb.jvarkit.util.vcf.VCFUtils)
VCFIterator(htsjdk.variant.vcf.VCFIterator)
Parameter(com.beust.jcommander.Parameter)
NullOuputStream(com.github.lindenb.jvarkit.io.NullOuputStream)
AnnPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.AnnPredictionParser)
ParametersDelegate(com.beust.jcommander.ParametersDelegate)
ArrayList(java.util.ArrayList)
HashSet(java.util.HashSet)
BiPredicate(java.util.function.BiPredicate)
ContigDictComparator(com.github.lindenb.jvarkit.util.samtools.ContigDictComparator)
RuntimeIOException(htsjdk.samtools.util.RuntimeIOException)
PedigreeParser(com.github.lindenb.jvarkit.pedigree.PedigreeParser)
IOUtils(com.github.lindenb.jvarkit.io.IOUtils)
Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher)
VepPredictionParser(com.github.lindenb.jvarkit.util.vcf.predictions.VepPredictionParser)
SortingCollection(htsjdk.samtools.util.SortingCollection)
VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine)
VCFHeaderLineType(htsjdk.variant.vcf.VCFHeaderLineType)
JexlGenotypePredicate(com.github.lindenb.jvarkit.util.vcf.JexlGenotypePredicate)
Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree)
SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary)
IOException(java.io.IOException)
EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator)
Comparator(java.util.Comparator)
Sample(com.github.lindenb.jvarkit.pedigree.Sample)
VariantContext(htsjdk.variant.variantcontext.VariantContext)
SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary)
EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator)
ContigDictComparator(com.github.lindenb.jvarkit.util.samtools.ContigDictComparator)
PedigreeParser(com.github.lindenb.jvarkit.pedigree.PedigreeParser)
NullOuputStream(com.github.lindenb.jvarkit.io.NullOuputStream)
VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine)
VCFHeader(htsjdk.variant.vcf.VCFHeader)
HashSet(java.util.HashSet)
PrintWriter(java.io.PrintWriter)
RuntimeIOException(htsjdk.samtools.util.RuntimeIOException)
Sample(com.github.lindenb.jvarkit.pedigree.Sample)
RuntimeIOException(htsjdk.samtools.util.RuntimeIOException)
IOException(java.io.IOException)
VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine)
EOFException(java.io.EOFException)
RuntimeIOException(htsjdk.samtools.util.RuntimeIOException)
IOException(java.io.IOException)
VCFEncoder(htsjdk.variant.vcf.VCFEncoder)
GeneExtractorFactory(com.github.lindenb.jvarkit.util.vcf.predictions.GeneExtractorFactory)
Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree)
VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)
Example 100 with CloseableIterator
use of htsjdk.samtools.util.CloseableIterator in project jvarkit by lindenb.
the class VCFCompareGT method doWork.
@Override
public int doWork(final List<String> arguments) {
if (arguments.isEmpty()) {
LOG.error("VCFs missing.");
return -1;
}
final List<String> vcfLabelList;
if (!StringUtil.isBlank(this.vcfLabelsStr)) {
vcfLabelList = Arrays.stream(this.vcfLabelsStr.split("[,]")).filter(S -> !StringUtil.isBlank(S)).collect(Collectors.toList());
if (new HashSet<>(vcfLabelList).size() != arguments.size()) {
LOG.error("bad number of labels in : " + this.vcfLabelsStr);
return -1;
}
if (vcfLabelList.stream().anyMatch(S -> !S.matches("[0-9A-Za-z]+"))) {
LOG.error("bad label in : " + this.vcfLabelsStr);
return -1;
}
} else {
vcfLabelList = new ArrayList<>();
for (int i = 0; i < arguments.size(); ++i) {
vcfLabelList.add(String.valueOf(i + 1));
}
}
VariantComparator varcmp = new VariantComparator();
SortingCollection<Variant> variants = null;
final Set<String> sampleNames = new LinkedHashSet<>();
try {
variants = SortingCollection.newInstance(Variant.class, new VariantCodec(), varcmp, writingSortingCollection.getMaxRecordsInRam(), writingSortingCollection.getTmpPaths());
variants.setDestructiveIteration(true);
final Set<VCFHeaderLine> metaData = new HashSet<VCFHeaderLine>();
metaData.add(new VCFHeaderLine(getClass().getSimpleName(), "version:" + getVersion() + " command:" + getProgramCommandLine()));
final BiFunction<String, Integer, String> createName = (SN, IDX) -> SN + "_" + vcfLabelList.get(IDX);
for (int i = 0; i < arguments.size(); ++i) {
final File vcfFile = new File(arguments.get(i));
LOG.info("Opening " + vcfFile);
final VCFReader vcfFileReader = VCFReaderFactory.makeDefault().open(vcfFile, false);
final CloseableIterator<VariantContext> iter = vcfFileReader.iterator();
final VCFHeader header = vcfFileReader.getHeader();
sampleNames.addAll(header.getSampleNamesInOrder());
metaData.add(new VCFHeaderLine(getClass().getSimpleName() + "_" + vcfLabelList.get(i), "File: " + vcfFile.getPath()));
long nLines = 0;
while (iter.hasNext()) {
final VariantContext var = iter.next();
if (nLines++ % 10000 == 0) {
LOG.info(vcfFile + " " + nLines);
}
if (!this.variantFilter.test(var))
continue;
if (!var.isVariant())
continue;
if (!var.hasGenotypes())
continue;
for (final Genotype genotype : var.getGenotypes()) {
final Variant rec = new Variant();
rec.file_index1 = i + 1;
rec.sampleName = genotype.getSampleName();
rec.chrom = var.getContig();
rec.start = var.getStart();
rec.end = var.getEnd();
rec.ref = var.getReference().getDisplayString().toUpperCase();
if (var.hasID()) {
rec.id = var.getID();
}
if (genotype.hasDP()) {
rec.dp = genotype.getDP();
}
if (genotype.hasGQ()) {
rec.gq = genotype.getGQ();
}
final List<Allele> alleles = genotype.getAlleles();
if (genotype.isNoCall() || !genotype.isAvailable() || alleles == null) {
if (!this.convertNoCallToHomRef)
continue;
rec.a1 = var.getReference().getDisplayString().toUpperCase();
rec.a2 = rec.a1;
} else if (alleles.size() == 1) {
rec.a1 = alleles.get(0).getDisplayString().toUpperCase();
rec.a2 = rec.a1;
} else if (alleles.size() == 2) {
rec.a1 = alleles.get(0).getDisplayString().toUpperCase();
rec.a2 = alleles.get(1).getDisplayString().toUpperCase();
if (rec.a1.compareTo(rec.a2) > 0) {
final String tmp = rec.a2;
rec.a2 = rec.a1;
rec.a1 = tmp;
}
} else {
continue;
}
variants.add(rec);
}
}
iter.close();
vcfFileReader.close();
}
variants.doneAdding();
LOG.info("Done Adding");
final String GenpotypeChangedKey = "GCH";
final String GenpotypeCreated = "GNW";
final String GenpotypeDiff = "GDF";
metaData.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.GENOTYPE_KEY));
metaData.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.DEPTH_KEY));
metaData.add(VCFStandardHeaderLines.getFormatLine(VCFConstants.GENOTYPE_QUALITY_KEY));
metaData.add(new VCFFormatHeaderLine(GenpotypeChangedKey, 1, VCFHeaderLineType.Integer, "Changed Genotype"));
metaData.add(new VCFFormatHeaderLine(GenpotypeCreated, 1, VCFHeaderLineType.Integer, "Genotype Created/Deleted"));
metaData.add(new VCFInfoHeaderLine(GenpotypeDiff, VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, "Samples with Genotype Difference"));
metaData.add(new VCFFilterHeaderLine("DISCORDANCE", "something has changed."));
final Set<String> newSampleNames = new TreeSet<>();
for (int i = 0; i < vcfLabelList.size(); ++i) {
for (final String sample : sampleNames) {
newSampleNames.add(createName.apply(sample, i));
}
}
final VCFHeader header = new VCFHeader(metaData, new ArrayList<>(newSampleNames));
final VariantContextWriter w = super.openVariantContextWriter(outputFile);
w.writeHeader(header);
final PosComparator posCompare = new PosComparator();
final EqualRangeIterator<Variant> iter = new EqualRangeIterator<>(variants.iterator(), posCompare);
while (iter.hasNext()) {
final List<Variant> row = iter.next();
/**
* this sample is not always the same
*/
final Set<String> samplesModified = new TreeSet<>();
/**
* the number of sample is different from vcflist.size()
*/
final Set<String> samplesCreates = new TreeSet<>();
final Map<String, List<Variant>> sample2variants = row.stream().collect(Collectors.groupingBy(T -> T.sampleName));
for (final String sn : sample2variants.keySet()) {
boolean all_hom_ref = true;
final List<Variant> sampleVariants = sample2variants.get(sn);
for (int x = 0; x < /*+1 non, besoin de tester hom_ref */
sampleVariants.size(); ++x) {
final Variant var1 = sampleVariants.get(x);
if (!var1.isHomRef())
all_hom_ref = false;
for (int y = x + 1; y < sampleVariants.size(); ++y) {
final Variant var2 = sampleVariants.get(y);
if (var1.a1.equals(var2.a1) && var1.a2.equals(var2.a2))
continue;
samplesModified.add(var1.sampleName);
}
}
if (sampleVariants.size() != arguments.size()) {
if (!convertNoCallToHomRef || (this.convertNoCallToHomRef && !all_hom_ref)) {
samplesCreates.add(sn);
}
}
}
final Variant first = row.get(0);
final Set<Allele> alleles = new HashSet<>();
alleles.add(Allele.create(first.ref, true));
for (final Variant var : row) {
alleles.add(Allele.create(var.a1, var.a1.equalsIgnoreCase(var.ref)));
alleles.add(Allele.create(var.a2, var.a2.equalsIgnoreCase(var.ref)));
}
final VariantContextBuilder b = new VariantContextBuilder(getClass().getName(), first.chrom, first.start, first.end, alleles);
// build genotypes
final List<Genotype> genotypes = new ArrayList<Genotype>();
for (final Variant var : row) {
// alleles for this genotype
final List<Allele> galleles = new ArrayList<Allele>();
galleles.add(Allele.create(var.a1, var.a1.equalsIgnoreCase(var.ref)));
galleles.add(Allele.create(var.a2, var.a2.equalsIgnoreCase(var.ref)));
final GenotypeBuilder gb = new GenotypeBuilder();
gb.DP(var.dp);
gb.alleles(galleles);
gb.name(createName.apply(var.sampleName, var.file_index1 - 1));
gb.GQ(var.gq);
gb.attribute(GenpotypeChangedKey, samplesModified.contains(var.sampleName) ? 1 : 0);
gb.attribute(GenpotypeCreated, samplesCreates.contains(var.sampleName) ? 1 : 0);
genotypes.add(gb.make());
}
b.genotypes(genotypes);
b.id(first.id);
if (!(samplesModified.isEmpty() && samplesCreates.isEmpty())) {
final Set<String> set2 = new TreeSet<String>(samplesModified);
set2.addAll(samplesCreates);
b.attribute(GenpotypeDiff, set2.toArray());
b.filter("DISCORDANCE");
} else {
b.passFilters();
}
if (!only_print_modified || !(samplesModified.isEmpty() && samplesCreates.isEmpty())) {
w.add(b.make());
}
}
iter.close();
w.close();
} catch (final Exception err) {
LOG.error(err);
return -1;
} finally {
if (variants != null)
try {
variants.cleanup();
} catch (Exception err) {
}
}
return 0;
}
Also used :
LinkedHashSet(java.util.LinkedHashSet)
Genotype(htsjdk.variant.variantcontext.Genotype)
VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine)
DataInputStream(java.io.DataInputStream)
CloseableIterator(htsjdk.samtools.util.CloseableIterator)
Allele(htsjdk.variant.variantcontext.Allele)
Arrays(java.util.Arrays)
JexlVariantPredicate(com.github.lindenb.jvarkit.util.vcf.JexlVariantPredicate)
Program(com.github.lindenb.jvarkit.util.jcommander.Program)
VCFStandardHeaderLines(htsjdk.variant.vcf.VCFStandardHeaderLines)
Parameter(com.beust.jcommander.Parameter)
VCFHeader(htsjdk.variant.vcf.VCFHeader)
BiFunction(java.util.function.BiFunction)
TreeSet(java.util.TreeSet)
ParametersDelegate(com.beust.jcommander.ParametersDelegate)
ArrayList(java.util.ArrayList)
HashSet(java.util.HashSet)
DataOutputStream(java.io.DataOutputStream)
StringUtil(htsjdk.samtools.util.StringUtil)
AbstractDataCodec(com.github.lindenb.jvarkit.util.picard.AbstractDataCodec)
Map(java.util.Map)
Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher)
VCFReaderFactory(com.github.lindenb.jvarkit.variant.vcf.VCFReaderFactory)
LinkedHashSet(java.util.LinkedHashSet)
VCFConstants(htsjdk.variant.vcf.VCFConstants)
SortingCollection(htsjdk.samtools.util.SortingCollection)
VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine)
VCFHeaderLineType(htsjdk.variant.vcf.VCFHeaderLineType)
GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder)
Predicate(java.util.function.Predicate)
Logger(com.github.lindenb.jvarkit.util.log.Logger)
VCFReader(htsjdk.variant.vcf.VCFReader)
Set(java.util.Set)
IOException(java.io.IOException)
Collectors(java.util.stream.Collectors)
File(java.io.File)
List(java.util.List)
VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter)
EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator)
VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine)
VariantContext(htsjdk.variant.variantcontext.VariantContext)
VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine)
Comparator(java.util.Comparator)
VCFHeaderLineCount(htsjdk.variant.vcf.VCFHeaderLineCount)
VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)
VCFHeaderLine(htsjdk.variant.vcf.VCFHeaderLine)
ArrayList(java.util.ArrayList)
VariantContext(htsjdk.variant.variantcontext.VariantContext)
EqualRangeIterator(com.github.lindenb.jvarkit.util.iterator.EqualRangeIterator)
VCFReader(htsjdk.variant.vcf.VCFReader)
TreeSet(java.util.TreeSet)
ArrayList(java.util.ArrayList)
List(java.util.List)
VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine)
VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter)
VCFHeader(htsjdk.variant.vcf.VCFHeader)
HashSet(java.util.HashSet)
LinkedHashSet(java.util.LinkedHashSet)
VCFFormatHeaderLine(htsjdk.variant.vcf.VCFFormatHeaderLine)
Genotype(htsjdk.variant.variantcontext.Genotype)
GenotypeBuilder(htsjdk.variant.variantcontext.GenotypeBuilder)
VCFInfoHeaderLine(htsjdk.variant.vcf.VCFInfoHeaderLine)
IOException(java.io.IOException)
Allele(htsjdk.variant.variantcontext.Allele)
VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)
File(java.io.File)
Aggregations
CloseableIterator (htsjdk.samtools.util.CloseableIterator)103
List (java.util.List)86
Logger (com.github.lindenb.jvarkit.util.log.Logger)85
Parameter (com.beust.jcommander.Parameter)82
Program (com.github.lindenb.jvarkit.util.jcommander.Program)78
ArrayList (java.util.ArrayList)73
Collectors (java.util.stream.Collectors)71
SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)69
Path (java.nio.file.Path)69
Launcher (com.github.lindenb.jvarkit.util.jcommander.Launcher)66
CloserUtil (htsjdk.samtools.util.CloserUtil)64
Set (java.util.Set)64
VCFHeader (htsjdk.variant.vcf.VCFHeader)59
VariantContext (htsjdk.variant.variantcontext.VariantContext)54
IOException (java.io.IOException)53
SequenceDictionaryUtils (com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils)51
SAMRecord (htsjdk.samtools.SAMRecord)51
VariantContextWriter (htsjdk.variant.variantcontext.writer.VariantContextWriter)51
IOUtils (com.github.lindenb.jvarkit.io.IOUtils)50
StringUtils (com.github.lindenb.jvarkit.lang.StringUtils)49
