Search in sources :

Example 31 with FastqRecord

use of htsjdk.samtools.fastq.FastqRecord in project gridss by PapenfussLab.

the class StructuralVariationCallBuilderTest method breakpoint_assembly_should_be_written.

@Test
@Ignore("Issue #17: Output breakpoint assembly sequences")
public void breakpoint_assembly_should_be_written() {
    ProcessingContext pc = getContext();
    String seq = "CATTAATCGCAAGAGCGGGTTGTATTCGcCGCCAAGTCAGCTGAAGCACCATTACCCGAtCAAAACATATCAGAAATGATTGACGTATCACAAGCCGGATTTTGTTTACAGCCTGTCTTATATCCTGAATAACGCACCGCCTATTCG";
    int anchor = 78;
    SAMRecord ass = AssemblyFactory.createAnchoredBreakend(getContext(), AES(), new SequentialIdGenerator("asm"), FWD, null, 6, 78, anchor, B(seq), B(40, seq.length()));
    SingleReadEvidence ae = incorporateRealignment(AES(), ass, ImmutableList.of(withQual(B(40, seq.length() - anchor), Read(2, anchor + 1, String.format("%dM", seq.length() - anchor)))[0]));
    StructuralVariationCallBuilder builder = new StructuralVariationCallBuilder(pc, (VariantContextDirectedEvidence) new IdsvVariantContextBuilder(getContext()) {

        {
            breakpoint(new BreakpointSummary(6, FWD, 78, 2, BWD, 79), "");
            phredScore(50);
        }
    }.make());
    builder.addEvidence(ae);
    VariantContextDirectedEvidence e = builder.make();
    FastqRecord fr = e.getBreakendAssemblyFastq();
    assertEquals(e.getEvidenceID() + "#" + Integer.toString(anchor), fr.getReadName());
    assertEquals(seq, fr.getReadString());
    assertEquals(fr.getReadName(), e.getAttribute(VcfSvConstants.BREAKPOINT_ID_KEY));
}
Also used : SAMRecord(htsjdk.samtools.SAMRecord) FastqRecord(htsjdk.samtools.fastq.FastqRecord) Ignore(org.junit.Ignore) Test(org.junit.Test)

Example 32 with FastqRecord

use of htsjdk.samtools.fastq.FastqRecord in project gridss by PapenfussLab.

the class SplitReadRealigner method processAlignmentRecord.

private void processAlignmentRecord(StreamingAligner aligner, SplitReadFastqExtractor recursiveExtractor, Map<String, SplitReadRealignmentInfo> realignments, SAMFileWriter writer) throws IOException {
    SAMRecord supp = aligner.getAlignment();
    String lookupkey = SplitReadIdentificationHelper.getOriginatingAlignmentUniqueName(supp);
    SplitReadRealignmentInfo info = realignments.get(lookupkey);
    if (supp.getSupplementaryAlignmentFlag() || supp.getNotPrimaryAlignmentFlag()) {
    // only consider the best mapping location reported by the aligner
    // log.debug(String.format("%s: ignoring supp alignment", supp.getReadName()));
    } else {
        assert (info.outstandingRealignments > 0);
        info.outstandingRealignments--;
        if (!supp.getReadUnmappedFlag()) {
            info.realignments.add(supp);
            List<FastqRecord> nestedRealignments = recursiveExtractor.extract(supp);
            for (FastqRecord fq : nestedRealignments) {
                aligner.asyncAlign(fq);
                info.outstandingRealignments++;
            // log.debug(String.format("%s: performing nested realignment. %d realignments now outstanding", info.originatingRecord.getReadName(), info.outstandingRealignments));
            }
        }
        // all splits identified
        if (info.outstandingRealignments == 0) {
            if (info.realignments.size() > 0) {
                SplitReadIdentificationHelper.convertToSplitRead(info.originatingRecord, info.realignments);
            }
            writer.addAlignment(info.originatingRecord);
            for (SAMRecord sar : info.realignments) {
                writer.addAlignment(sar);
            }
            // log.debug(String.format("%s: %d supp alignments found.", info.originatingRecord.getReadName(), info.realignments.size()));
            realignments.remove(lookupkey);
        } else {
        // log.debug(String.format("%s: %d outstanding alignments", info.originatingRecord.getReadName(), info.outstandingRealignments));
        }
    }
}
Also used : SAMRecord(htsjdk.samtools.SAMRecord) FastqRecord(htsjdk.samtools.fastq.FastqRecord)

Example 33 with FastqRecord

use of htsjdk.samtools.fastq.FastqRecord in project gridss by PapenfussLab.

the class SplitReadIdentificationHelper method getSplitReadRealignments.

/**
 * Extract the unaligned portions of the read requiring realignment to identify split reads
 *
 * @param r alignment record
 * @param recordIsPartialAlignment true if the record is the result of aligning FastqRecords
 * from a previous call to getSplitReadRealignments()
 * @return bases requiring alignment to identify split reads
 */
public static List<FastqRecord> getSplitReadRealignments(SAMRecord r, boolean recordIsPartialAlignment, EvidenceIdentifierGenerator eidgen) {
    int startClipLength = SAMRecordUtil.getStartSoftClipLength(r);
    int endClipLength = SAMRecordUtil.getEndSoftClipLength(r);
    if (startClipLength + endClipLength == 0 || r.getReadUnmappedFlag()) {
        return Collections.emptyList();
    }
    int offset;
    String name;
    if (recordIsPartialAlignment) {
        offset = getRealignmentFirstAlignedBaseReadOffset(r);
        name = getOriginatingAlignmentUniqueName(r);
    } else {
        offset = 0;
        name = eidgen.getAlignmentUniqueName(r);
    }
    List<FastqRecord> list = new ArrayList<>(2);
    if (startClipLength > 0) {
        byte[] startbases = Arrays.copyOfRange(r.getReadBases(), 0, startClipLength);
        byte[] startqual = Arrays.copyOfRange(r.getBaseQualities(), 0, startClipLength);
        int startOffset = offset;
        if (r.getReadNegativeStrandFlag()) {
            SequenceUtil.reverseComplement(startbases);
            ArrayUtils.reverse(startqual);
            // 5432
            // SS
            startOffset = startOffset + r.getReadLength() - startClipLength;
        }
        FastqRecord start = new FastqRecord(getSplitAlignmentFastqName(name, startOffset), new String(startbases, StandardCharsets.US_ASCII), "", SAMUtils.phredToFastq(startqual));
        list.add(start);
    }
    if (endClipLength > 0) {
        byte[] endbases = Arrays.copyOfRange(r.getReadBases(), r.getReadLength() - endClipLength, r.getReadLength());
        byte[] endqual = Arrays.copyOfRange(r.getBaseQualities(), r.getReadLength() - endClipLength, r.getReadLength());
        int endOffset = offset + r.getReadLength() - endClipLength;
        if (r.getReadNegativeStrandFlag()) {
            SequenceUtil.reverseComplement(endbases);
            ArrayUtils.reverse(endqual);
            // 5432
            // SS
            endOffset = offset;
        }
        FastqRecord end = new FastqRecord(getSplitAlignmentFastqName(name, endOffset), new String(endbases, StandardCharsets.US_ASCII), "", SAMUtils.phredToFastq(endqual));
        list.add(end);
    }
    return list;
}
Also used : FastqRecord(htsjdk.samtools.fastq.FastqRecord) ArrayList(java.util.ArrayList)

Example 34 with FastqRecord

use of htsjdk.samtools.fastq.FastqRecord in project jvarkit by lindenb.

the class FastqRevComp method run.

private void run(FastqReader r, PrintStream out) {
    String s;
    long nRec = 0L;
    r.setValidationStringency(ValidationStringency.LENIENT);
    while (r.hasNext()) {
        if (++nRec % 1E6 == 0) {
            LOG.info("N-Reads:" + nRec);
        }
        FastqRecord fastq = r.next();
        out.print(FastqConstants.SEQUENCE_HEADER);
        out.println(fastq.getReadName());
        s = fastq.getReadString();
        if (// interleaced
        (this.only_R2 && nRec % 2 == 1) || (this.only_R1 && nRec % 2 == 0)) {
            out.print(s);
        } else {
            for (int i = s.length() - 1; i >= 0; i--) {
                out.print(AcidNucleics.complement(s.charAt(i)));
            }
        }
        out.println();
        out.print(FastqConstants.QUALITY_HEADER);
        s = fastq.getBaseQualityHeader();
        if (s != null)
            out.print(s);
        out.println();
        s = fastq.getBaseQualityString();
        if (// interleaced
        (this.only_R2 && nRec % 2 == 1) || (this.only_R1 && nRec % 2 == 0)) {
            out.print(s);
        } else {
            for (int i = s.length() - 1; i >= 0; i--) {
                out.print(s.charAt(i));
            }
        }
        out.println();
        if (out.checkError())
            break;
    }
    out.flush();
    LOG.info("Done. N-Reads:" + nRec);
}
Also used : FastqRecord(htsjdk.samtools.fastq.FastqRecord)

Example 35 with FastqRecord

use of htsjdk.samtools.fastq.FastqRecord in project jvarkit by lindenb.

the class FastqToFasta method run.

private void run(FastqReader r, PrintStream out) {
    int wsp = 0;
    long nRec = 0L;
    r.setValidationStringency(ValidationStringency.LENIENT);
    while (r.hasNext()) {
        if (++nRec % 1E6 == 0) {
            LOG.info("N-Reads:" + nRec);
        }
        FastqRecord fastq = r.next();
        out.print(">");
        if (!trim_after_space || (wsp = fastq.getReadName().indexOf(' ')) == -1) {
            out.println(fastq.getReadName());
        } else {
            out.println(fastq.getReadName().substring(0, wsp));
        }
        int readLen = fastq.getReadString().length();
        int i = 0;
        while (i < readLen) {
            int end = Math.min(i + fastaLineLen, readLen);
            out.println(fastq.getReadString().substring(i, end));
            i = end;
        }
        if (out.checkError())
            break;
    }
    out.flush();
    LOG.info("Done. N-Reads:" + nRec);
}
Also used : FastqRecord(htsjdk.samtools.fastq.FastqRecord)

Aggregations

FastqRecord (htsjdk.samtools.fastq.FastqRecord)40 SAMRecord (htsjdk.samtools.SAMRecord)16 Test (org.junit.Test)12 FastqReader (htsjdk.samtools.fastq.FastqReader)8 File (java.io.File)6 SAMFileHeader (htsjdk.samtools.SAMFileHeader)4 FastqReader (com.github.lindenb.jvarkit.util.picard.FastqReader)3 FourLinesFastqReader (com.github.lindenb.jvarkit.util.picard.FourLinesFastqReader)3 SAMFileWriter (htsjdk.samtools.SAMFileWriter)3 IOException (java.io.IOException)3 Test (org.testng.annotations.Test)3 BufferedReferenceSequenceFile (au.edu.wehi.idsv.picard.BufferedReferenceSequenceFile)2 ChimericAlignment (au.edu.wehi.idsv.sam.ChimericAlignment)2 SAMSequenceDictionaryProgress (com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress)2 Cigar (htsjdk.samtools.Cigar)2 CigarElement (htsjdk.samtools.CigarElement)2 SAMFileWriterFactory (htsjdk.samtools.SAMFileWriterFactory)2 SAMRecordIterator (htsjdk.samtools.SAMRecordIterator)2 CommandLineProgramTest (org.broadinstitute.hellbender.CommandLineProgramTest)2 ProgressLogger (org.broadinstitute.hellbender.utils.runtime.ProgressLogger)2