use of com.github.lindenb.jvarkit.util.log.Logger in project jvarkit by lindenb.
the class Biostar145820 method doWork.
@Override
public int doWork(final List<String> args) {
SamReader samReader = null;
SAMRecordIterator iter = null;
SAMFileWriter samWriter = null;
Random random = new Random();
CloseableIterator<RandSamRecord> iter2 = null;
try {
final String input = oneFileOrNull(args);
samReader = super.openSamReader(input);
final SAMFileHeader header = samReader.getFileHeader().clone();
header.setSortOrder(SortOrder.unsorted);
header.addComment("Processed with " + getProgramName() + " : " + getProgramCommandLine());
final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(samReader.getFileHeader()).logger(LOG);
iter = samReader.iterator();
final SortingCollection<RandSamRecord> sorter = SortingCollection.newInstance(RandSamRecord.class, new RandSamRecordCodec(header), new RandSamRecordComparator(), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
sorter.setDestructiveIteration(true);
while (iter.hasNext()) {
final SAMRecord rec = progress.watch(iter.next());
if (this.filter.filterOut(rec)) {
continue;
}
final RandSamRecord r = new RandSamRecord();
r.rand_index = random.nextInt();
r.samRecord = progress.watch(rec);
sorter.add(r);
}
iter.close();
iter = null;
sorter.doneAdding();
iter2 = sorter.iterator();
samWriter = writingBamArgs.openSAMFileWriter(outputFile, header, true);
final SAMFileWriter finalw = samWriter;
iter2.stream().limit(this.count < 0L ? Long.MAX_VALUE - 1 : this.count).map(R -> R.samRecord).forEach(R -> finalw.addAlignment(R));
iter2.close();
iter2 = null;
sorter.cleanup();
progress.finish();
} catch (final Exception e) {
LOG.error(e);
return -1;
} finally {
CloserUtil.close(iter);
CloserUtil.close(iter2);
CloserUtil.close(samReader);
CloserUtil.close(samWriter);
}
return 0;
}
use of com.github.lindenb.jvarkit.util.log.Logger in project jvarkit by lindenb.
the class VcfLoopOverGenes method doWork.
@SuppressWarnings("resource")
@Override
public int doWork(final List<String> args) {
PrintWriter pw = null;
VCFFileReader vcfFileReader = null;
CloseableIterator<VariantContext> iter = null;
CloseableIterator<GeneLoc> iter2 = null;
BufferedReader br = null;
ArchiveFactory archive = null;
try {
final File vcf = new File(oneAndOnlyOneFile(args));
vcfFileReader = new VCFFileReader(vcf, (this.geneFile != null || !StringUtil.isBlank(this.regionStr)));
this.dictionary = vcfFileReader.getFileHeader().getSequenceDictionary();
if (this.dictionary == null) {
throw new JvarkitException.VcfDictionaryMissing(vcf);
}
final VcfTools tools = new VcfTools(vcfFileReader.getFileHeader());
if (!this.prefix.isEmpty() && !this.prefix.endsWith(".")) {
this.prefix += ".";
}
if (this.geneFile == null) {
final SortingCollection<GeneLoc> sortingCollection = SortingCollection.newInstance(GeneLoc.class, new GeneLocCodec(), (A, B) -> A.compareTo(B), this.writingSortingCollection.getMaxRecordsInRam(), this.writingSortingCollection.getTmpPaths());
sortingCollection.setDestructiveIteration(true);
if (StringUtil.isBlank(this.regionStr)) {
iter = vcfFileReader.iterator();
} else {
final IntervalParser parser = new IntervalParser(this.dictionary);
parser.setContigNameIsWholeContig(true);
final Interval interval = parser.parse(this.regionStr);
if (interval == null) {
LOG.error("Cannot parse interval " + this.regionStr);
return -1;
}
iter = vcfFileReader.query(interval.getContig(), interval.getStart(), interval.getEnd());
}
final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(vcfFileReader.getFileHeader()).logger(LOG);
if (this.splitMethod.equals(SplitMethod.Annotations)) {
while (iter.hasNext()) {
final VariantContext ctx = progress.watch(iter.next());
for (final AnnPredictionParser.AnnPrediction pred : tools.getAnnPredictionParser().getPredictions(ctx)) {
if (this.snpEffNoIntergenic && pred.isIntergenicRegion()) {
continue;
}
if (!StringUtil.isBlank(pred.getGeneName())) {
sortingCollection.add(create(ctx, pred.getGeneName(), SourceType.ANN_GeneName));
}
if (!StringUtil.isBlank(pred.getGeneId())) {
sortingCollection.add(create(ctx, pred.getGeneId(), SourceType.ANN_GeneID));
}
if (!StringUtil.isBlank(pred.getFeatureId())) {
sortingCollection.add(create(ctx, pred.getFeatureId(), SourceType.ANN_FeatureID));
}
}
for (final VepPredictionParser.VepPrediction pred : tools.getVepPredictionParser().getPredictions(ctx)) {
if (!StringUtil.isBlank(pred.getGene())) {
sortingCollection.add(create(ctx, pred.getGene(), SourceType.VEP_Gene));
}
if (!StringUtil.isBlank(pred.getFeature())) {
sortingCollection.add(create(ctx, pred.getFeature(), SourceType.VEP_Feature));
}
if (!StringUtil.isBlank(pred.getSymbol())) {
sortingCollection.add(create(ctx, pred.getSymbol(), SourceType.VEP_Symbol));
}
if (!StringUtil.isBlank(pred.getHgncId())) {
sortingCollection.add(create(ctx, pred.getHgncId(), SourceType.VEP_HgncId));
}
}
}
} else /**
* split VCF per sliding window of variants
*/
if (this.splitMethod.equals(SplitMethod.VariantSlidingWindow)) {
if (this.variantsWinCount < 1) {
LOG.error("Bad value for variantsWinCount");
return -1;
}
if (this.variantsWinShift < 1 || this.variantsWinShift > this.variantsWinCount) {
LOG.error("Bad value for variantsWinShift");
return -1;
}
final List<VariantContext> buffer = new ArrayList<>(this.variantsWinCount);
/**
* routine to dump buffer into sorting collection
*/
final Runnable dumpBuffer = () -> {
if (buffer.isEmpty())
return;
final String contig = buffer.get(0).getContig();
final int chromStart = buffer.stream().mapToInt(CTX -> CTX.getStart()).min().getAsInt();
// use last of start too
final int chromEnd0 = buffer.stream().mapToInt(CTX -> CTX.getStart()).max().getAsInt();
// final int chromEnd1 = buffer.stream().mapToInt(CTX->CTX.getEnd()).max().getAsInt();
final String identifier = contig + "_" + String.format(NUM_FORMAT, chromStart) + "_" + String.format(NUM_FORMAT, chromEnd0);
for (final VariantContext ctx : buffer) {
sortingCollection.add(create(ctx, identifier, SourceType.SlidingVariants));
}
};
while (iter.hasNext()) {
VariantContext ctx = progress.watch(iter.next());
/* reduce the memory footprint for this context */
ctx = new VariantContextBuilder(ctx).genotypes(Collections.emptyList()).unfiltered().rmAttributes(new ArrayList<>(ctx.getAttributes().keySet())).make();
if (!buffer.isEmpty() && !buffer.get(0).getContig().equals(ctx.getContig())) {
dumpBuffer.run();
buffer.clear();
}
buffer.add(ctx);
if (buffer.size() >= this.variantsWinCount) {
dumpBuffer.run();
final int fromIndex = Math.min(this.variantsWinShift, buffer.size());
buffer.subList(0, fromIndex).clear();
}
}
dumpBuffer.run();
buffer.clear();
} else if (this.splitMethod.equals(SplitMethod.ContigSlidingWindow)) {
if (this.contigWinLength < 1) {
LOG.error("Bad value for contigWinCount");
return -1;
}
if (this.contigWinShift < 1 || this.contigWinShift > this.contigWinLength) {
LOG.error("Bad value for contigWinShift");
return -1;
}
while (iter.hasNext()) {
VariantContext ctx = progress.watch(iter.next());
/* reduce the memory footprint for this context */
ctx = new VariantContextBuilder(ctx).genotypes(Collections.emptyList()).unfiltered().rmAttributes(new ArrayList<>(ctx.getAttributes().keySet())).make();
int start = 0;
while (start <= ctx.getStart()) {
if (start + this.contigWinLength >= ctx.getStart()) {
final int chromStart = start;
final int chromEnd0 = start + this.contigWinLength;
final String identifier = ctx.getContig() + "_" + String.format(NUM_FORMAT, chromStart) + "_" + String.format(NUM_FORMAT, chromEnd0);
sortingCollection.add(create(ctx, identifier, SourceType.SlidingContig));
}
start += this.contigWinShift;
}
}
} else {
throw new IllegalStateException("No such method: " + this.splitMethod);
}
sortingCollection.doneAdding();
progress.finish();
iter.close();
iter = null;
pw = super.openFileOrStdoutAsPrintWriter(this.outputFile);
iter2 = sortingCollection.iterator();
final EqualRangeIterator<GeneLoc> eqiter = new EqualRangeIterator<>(iter2, this.compareGeneName);
int geneIdentifierId = 0;
while (eqiter.hasNext()) {
final List<GeneLoc> gene = eqiter.next();
pw.print(gene.get(0).contig);
pw.print('\t');
// -1 for BED
pw.print(gene.stream().mapToInt(G -> G.start).min().getAsInt() - 1);
pw.print('\t');
pw.print(gene.stream().mapToInt(G -> G.end).max().getAsInt());
pw.print('\t');
pw.print(this.prefix + String.format("%09d", ++geneIdentifierId));
pw.print('\t');
pw.print(gene.get(0).geneName);
pw.print('\t');
pw.print(gene.get(0).sourceType);
pw.print('\t');
pw.print(gene.size());
pw.println();
}
pw.flush();
pw.close();
pw = null;
eqiter.close();
iter2.close();
iter2 = null;
sortingCollection.cleanup();
} else {
if (this.nJobs < 1) {
this.nJobs = Math.max(1, Runtime.getRuntime().availableProcessors());
LOG.info("setting njobs to " + this.nJobs);
}
final ExecutorService executorService;
final List<Future<Integer>> futureResults;
if (this.nJobs > 1) {
executorService = new ThreadPoolExecutor(this.nJobs, this.nJobs, 0L, TimeUnit.MILLISECONDS, new LinkedBlockingQueue<Runnable>());
futureResults = new ArrayList<>();
} else {
executorService = null;
futureResults = Collections.emptyList();
}
if (this.outputFile == null) {
LOG.error("When scanning a VCF with " + this.geneFile + ". Output file must be defined");
}
if (!this.exec.isEmpty()) {
if (this.outputFile.getName().endsWith(".zip")) {
LOG.error("Cannot execute " + this.exec + " when saving to a zip.");
return -1;
}
}
archive = ArchiveFactory.open(this.outputFile);
PrintWriter manifest = this.deleteAfterCommand && !this.exec.isEmpty() ? // all files will be deleted, no manifest needed
new PrintWriter(new NullOuputStream()) : archive.openWriter(this.prefix + "manifest.txt");
br = IOUtils.openFileForBufferedReading(this.geneFile);
final BedLineCodec bedCodec = new BedLineCodec();
for (; ; ) {
if (!futureResults.isEmpty()) {
int i = 0;
while (i < futureResults.size()) {
final Future<Integer> r = futureResults.get(i);
if (r.isCancelled()) {
LOG.error("Task was canceled. Break.");
return -1;
} else if (r.isDone()) {
futureResults.remove(i);
int rez = r.get();
if (rez != 0) {
LOG.error("Task Failed (" + rez + "). Break");
}
} else {
i++;
}
}
}
final String line = br.readLine();
if (line == null)
break;
if (line.startsWith("#") || line.isEmpty())
continue;
final BedLine bedLine = bedCodec.decode(line);
if (bedLine == null)
continue;
// ID
final String geneIdentifier = bedLine.get(3);
// name
final String geneName = bedLine.get(4);
final SourceType sourceType = SourceType.valueOf(bedLine.get(5));
final String filename = geneIdentifier;
final String outputVcfName = (filename.startsWith(this.prefix) ? "" : this.prefix) + filename + ".vcf" + (this.compress ? ".gz" : "");
LOG.info(bedLine.getContig() + ":" + bedLine.getStart() + "-" + bedLine.getEnd() + " length :" + (bedLine.getEnd() - bedLine.getStart()));
if (bedLine.getEnd() - bedLine.getStart() > 1E6) {
LOG.warn("That's a large region ! " + bedLine);
}
OutputStream vcfOutputStream = null;
VariantContextWriter vw = null;
int countVariants = 0;
final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(vcfFileReader.getFileHeader()).logger(LOG).prefix(geneName + " " + bedLine.getContig() + ":" + bedLine.getStart() + "-" + bedLine.getEnd());
iter = vcfFileReader.query(bedLine.getContig(), bedLine.getStart(), bedLine.getEnd());
while (iter.hasNext()) {
VariantContext ctx = progress.watch(iter.next());
switch(sourceType) {
case SlidingVariants:
{
// nothing
break;
}
case SlidingContig:
{
// nothing
break;
}
case ANN_GeneName:
case ANN_FeatureID:
case ANN_GeneID:
{
final List<String> preds = new ArrayList<>();
for (final AnnPredictionParser.AnnPrediction pred : tools.getAnnPredictionParser().getPredictions(ctx)) {
final String predictionIdentifier;
switch(sourceType) {
case ANN_GeneName:
predictionIdentifier = pred.getGeneName();
break;
case ANN_FeatureID:
predictionIdentifier = pred.getFeatureId();
break;
case ANN_GeneID:
predictionIdentifier = pred.getGeneId();
break;
default:
throw new IllegalStateException(bedLine.toString());
}
if (StringUtil.isBlank(predictionIdentifier))
continue;
if (!geneName.equals(predictionIdentifier))
continue;
preds.add(pred.getOriginalAttributeAsString());
}
if (preds.isEmpty()) {
ctx = null;
} else {
ctx = new VariantContextBuilder(ctx).rmAttribute(tools.getAnnPredictionParser().getTag()).attribute(tools.getAnnPredictionParser().getTag(), preds).make();
}
break;
}
case VEP_Gene:
case VEP_Feature:
case VEP_Symbol:
case VEP_HgncId:
{
final List<String> preds = new ArrayList<>();
for (final VepPredictionParser.VepPrediction pred : tools.getVepPredictions(ctx)) {
final String predictionIdentifier;
switch(sourceType) {
case VEP_Gene:
predictionIdentifier = pred.getGene();
break;
case VEP_Feature:
predictionIdentifier = pred.getFeature();
break;
case VEP_Symbol:
predictionIdentifier = pred.getSymbol();
break;
case VEP_HgncId:
predictionIdentifier = pred.getHgncId();
break;
default:
throw new IllegalStateException(bedLine.toString());
}
if (StringUtil.isBlank(predictionIdentifier))
continue;
if (!geneName.equals(predictionIdentifier))
continue;
preds.add(pred.getOriginalAttributeAsString());
}
if (preds.isEmpty()) {
ctx = null;
} else {
ctx = new VariantContextBuilder(ctx).rmAttribute(tools.getVepPredictionParser().getTag()).attribute(tools.getVepPredictionParser().getTag(), preds).make();
}
break;
}
default:
throw new IllegalStateException(bedLine.toString());
}
if (ctx == null)
continue;
if (vcfOutputStream == null) {
LOG.info(filename);
manifest.println(outputVcfName);
final VCFHeader header = new VCFHeader(vcfFileReader.getFileHeader());
header.addMetaDataLine(new VCFHeaderLine(VCF_HEADER_SPLITKEY, filename));
vcfOutputStream = archive.openOuputStream(outputVcfName);
vw = VCFUtils.createVariantContextWriterToOutputStream(vcfOutputStream);
vw.writeHeader(header);
}
countVariants++;
vw.add(ctx);
if (countVariants % 1000 == 0) {
LOG.info("Loading : " + geneIdentifier + " N=" + countVariants);
}
}
progress.finish();
LOG.info(geneIdentifier + " N=" + countVariants);
if (vcfOutputStream != null) {
vw.close();
vcfOutputStream.flush();
vcfOutputStream.close();
vw = null;
if (!this.exec.isEmpty()) {
final Callable<Integer> callable = () -> {
final File vcfOutFile = new File(this.outputFile, outputVcfName);
IOUtil.assertFileIsReadable(vcfOutFile);
final String vcfPath = vcfOutFile.getPath();
final StringTokenizer st = new StringTokenizer(this.exec);
final List<String> command = new ArrayList<>(1 + st.countTokens());
while (st.hasMoreTokens()) {
String token = st.nextToken().replaceAll("__PREFIX__", this.prefix).replaceAll("__CONTIG__", bedLine.getContig()).replaceAll("__CHROM__", bedLine.getContig()).replaceAll("__ID__", geneIdentifier).replaceAll("__NAME__", geneName).replaceAll("__START__", String.valueOf(bedLine.getStart())).replaceAll("__END__", String.valueOf(bedLine.getEnd())).replaceAll("__SOURCE__", sourceType.name()).replaceAll("__VCF__", vcfPath);
command.add(token);
}
LOG.info(command.stream().map(S -> "'" + S + "'").collect(Collectors.joining(" ")));
final ProcessBuilder pb = new ProcessBuilder(command);
pb.redirectErrorStream(true);
final Process p = pb.start();
final Thread stdoutThread = new Thread(() -> {
try {
InputStream in = p.getInputStream();
IOUtils.copyTo(in, stdout());
} catch (Exception err) {
LOG.error(err);
}
});
stdoutThread.start();
int exitValue = p.waitFor();
if (exitValue != 0) {
LOG.error("Command failed (" + exitValue + "):" + String.join(" ", command));
return -1;
} else {
if (deleteAfterCommand) {
if (!vcfOutFile.delete()) {
LOG.warn("Cannot delete " + vcfOutFile);
}
}
return 0;
}
};
if (executorService != null) {
final Future<Integer> rez = executorService.submit(callable);
futureResults.add(rez);
} else {
final int ret = callable.call();
if (ret != 0) {
LOG.error("Error with process (" + ret + ")");
return ret;
}
}
}
} else {
manifest.println("#" + filename);
LOG.warn("No Variant Found for " + line);
}
iter.close();
}
;
if (executorService != null) {
LOG.info("shutdown");
executorService.shutdown();
executorService.awaitTermination(365, TimeUnit.DAYS);
}
br.close();
br = null;
manifest.close();
archive.close();
archive = null;
LOG.info("Done");
}
vcfFileReader.close();
vcfFileReader = null;
return 0;
} catch (Exception e) {
LOG.error(e);
return -1;
} finally {
{
CloserUtil.close(iter2);
CloserUtil.close(iter);
CloserUtil.close(pw);
CloserUtil.close(vcfFileReader);
CloserUtil.close(br);
CloserUtil.close(archive);
}
}
}
use of com.github.lindenb.jvarkit.util.log.Logger in project jvarkit by lindenb.
the class VcfBurdenGoEnrichment method doWork.
@Override
public int doWork(final List<String> args) {
if (StringUtil.isBlank(this.readingGo.goUri)) {
LOG.error("Undefined GOs uri.");
return -1;
}
if (this.geneFile == null || !this.geneFile.exists()) {
LOG.error("Undefined gene file option.");
return -1;
}
try {
final GoTree gotree = this.readingGo.createParser().setIgnoreDbXRef(true).parse(this.readingGo.goUri);
List<GoTree.Term> terms = new ArrayList<>(gotree.getTerms());
final Map<GoTree.Term, Node> term2node = new HashMap<>();
// build the node TREE
while (!terms.isEmpty()) {
int i = 0;
while (i < terms.size()) {
final GoTree.Term t = terms.get(i);
if (!t.hasRelations()) {
term2node.put(t, new Node(t));
terms.remove(i);
} else if (t.getRelations().stream().allMatch(L -> term2node.containsKey(L.getTo()))) {
final Node n = new Node(t);
n.parents.addAll(t.getRelations().stream().map(L -> term2node.get(L.getTo())).collect(Collectors.toSet()));
term2node.put(t, n);
terms.remove(i);
} else {
i++;
}
}
}
terms = null;
final Set<String> unknownAcn = new HashSet<>();
final Map<String, Set<Node>> gene2node = new HashMap<>();
final BufferedReader r = IOUtils.openFileForBufferedReading(this.geneFile);
String line;
while ((line = r.readLine()) != null) {
if (line.isEmpty() || line.startsWith("#"))
continue;
final int t = line.indexOf('\t');
if (t == -1) {
r.close();
LOG.error("tab missing in " + line + " of " + this.geneFile);
return -1;
}
final String gene = line.substring(0, t).trim();
if (StringUtil.isBlank(gene)) {
r.close();
LOG.error("Emtpy gene in " + line);
return -1;
}
// using getTermByName because found sysnonym in GOA
final String termAcn = line.substring(t + 1).trim();
if (unknownAcn.contains(termAcn))
continue;
final GoTree.Term term = gotree.getTermByName(termAcn);
if (term == null && !unknownAcn.contains(termAcn)) {
unknownAcn.add(termAcn);
LOG.warning("Don't know this GO term in " + line + " of " + this.geneFile + ". Could be obsolete, synonym, go specific division. Skipping.");
continue;
}
final Node node = term2node.get(term);
if (node == null) {
r.close();
LOG.error("Don't know this node in " + line + " of " + this.geneFile);
return -1;
}
Set<Node> nodes = gene2node.get(gene);
if (nodes == null) {
nodes = new HashSet<>();
gene2node.put(gene, nodes);
}
node.numGenes++;
nodes.add(node);
}
;
// clean up
unknownAcn.clear();
r.close();
final VcfIterator iter = openVcfIterator(oneFileOrNull(args));
final VCFHeader header = iter.getHeader();
final VepPredictionParser vepParser = new VepPredictionParserFactory(header).get();
final AnnPredictionParser annParser = new AnnPredictionParserFactory(header).get();
final Set<Pedigree.Person> persons;
if (this.pedFile != null) {
final Pedigree pedigree = Pedigree.newParser().parse(this.pedFile);
persons = new Pedigree.CaseControlExtractor().extract(header, pedigree);
} else {
persons = new Pedigree.CaseControlExtractor().extract(header);
}
final Set<Pedigree.Person> affected = persons.stream().filter(P -> P.isAffected()).collect(Collectors.toSet());
final Set<Pedigree.Person> unaffected = persons.stream().filter(P -> P.isUnaffected()).collect(Collectors.toSet());
if (affected.isEmpty()) {
LOG.error("No Affected individual");
return -1;
}
if (unaffected.isEmpty()) {
LOG.error("No unaffected individual");
return -1;
}
final List<String> lookColumns = Arrays.asList("CCDS", "Feature", "ENSP", "Gene", "HGNC", "HGNC_ID", "SYMBOL", "RefSeq");
final Predicate<Genotype> isWildGenotype = G -> {
if (G == null)
return false;
return G.isHomRef();
};
final Predicate<Genotype> isAltGenotype = G -> {
if (G == null)
return false;
return G.isCalled() && !G.isHomRef();
};
final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(header).logger(LOG);
while (iter.hasNext()) {
final VariantContext ctx = progress.watch(iter.next());
if (!this.variantFilter.test(ctx))
continue;
final Set<String> genes = new HashSet<>();
for (final String predStr : ctx.getAttributeAsList(vepParser.getTag()).stream().map(O -> String.class.cast(O)).collect(Collectors.toList())) {
final VepPredictionParser.VepPrediction pred = vepParser.parseOnePrediction(ctx, predStr);
for (final String col : lookColumns) {
final String token = pred.getByCol(col);
if (!StringUtil.isBlank(token)) {
genes.add(token);
}
}
}
for (final String predStr : ctx.getAttributeAsList(annParser.getTag()).stream().map(O -> String.class.cast(O)).collect(Collectors.toList())) {
final AnnPredictionParser.AnnPrediction pred = annParser.parseOnePrediction(predStr);
final String token = pred.getGeneName();
if (!StringUtil.isBlank(token)) {
genes.add(token);
}
}
if (genes.isEmpty())
continue;
final Set<Node> nodes = genes.stream().filter(G -> gene2node.containsKey(G)).flatMap(G -> gene2node.get(G).stream()).collect(Collectors.toSet());
if (nodes.isEmpty())
continue;
final long unaffected_alt = unaffected.stream().map(P -> ctx.getGenotype(P.getId())).filter(G -> this.genotypeFilter.test(ctx, G)).filter(isAltGenotype).count();
final long affected_alt = affected.stream().map(P -> ctx.getGenotype(P.getId())).filter(G -> this.genotypeFilter.test(ctx, G)).filter(isAltGenotype).count();
/* no informative */
if (unaffected_alt + affected_alt == 0L) {
continue;
}
final long affected_ref = affected.stream().map(P -> ctx.getGenotype(P.getId())).filter(G -> this.genotypeFilter.test(ctx, G)).filter(isWildGenotype).count();
final long unaffected_ref = unaffected.stream().map(P -> ctx.getGenotype(P.getId())).filter(G -> this.genotypeFilter.test(ctx, G)).filter(isWildGenotype).count();
nodes.stream().forEach(N -> N.resetVisitedFlag());
nodes.stream().forEach(N -> N.visit(unaffected_ref, unaffected_alt, affected_ref, affected_alt));
}
iter.close();
progress.finish();
LOG.info("Calculating Fisher and dumping.. please wait");
final PrintWriter pw = super.openFileOrStdoutAsPrintWriter(this.outputFile);
pw.println("#go_term\tfisher\tname\tgo_term_depth\tcount_genes_in_this_node" + "\tunaffected_ref_gt" + "\tunaffected_alt_gt" + "\taffected_ref_gt" + "\taffected_alt_gt");
term2node.values().stream().filter(N -> this.show_never_seeen_term || N.sum() > 0L).sorted((n1, n2) -> Double.compare(n1.fisher(), n2.fisher())).forEach(N -> {
pw.print(N.goTerm.getAcn());
pw.print('\t');
pw.print(N.fisher());
pw.print("\t");
pw.print(N.goTerm.getName().replaceAll("[ \',\\-]+", "_"));
pw.print("\t");
pw.print(N.goTerm.getMinDepth());
pw.print('\t');
pw.print(N.numGenes);
pw.print('\t');
pw.print(N.unaffected_ref);
pw.print('\t');
pw.print(N.unaffected_alt);
pw.print('\t');
pw.print(N.affected_ref);
pw.print('\t');
pw.print(N.affected_alt);
pw.println();
});
pw.flush();
pw.close();
return 0;
} catch (final Exception err) {
LOG.error(err);
return -1;
}
}
Aggregations