Search in sources :

Example 6 with ChimericAlignment

use of au.edu.wehi.idsv.sam.ChimericAlignment in project gridss by PapenfussLab.

the class ExtractSVReadsTest method regression_should_extract_split_read_alignments_as_group_2.

// @Test
@Category(Hg38Tests.class)
public void regression_should_extract_split_read_alignments_as_group_2() throws IOException {
    File ref = Hg38Tests.findHg38Reference();
    ReferenceLookup lookup = new SynchronousReferenceLookupAdapter(new IndexedFastaSequenceFile(ref));
    Files.copy(new File("src/test/resources/sa.split/test2.sam"), input);
    ExtractSVReads extract = new ExtractSVReads();
    // new ProcessingContext(getFSContext(), ref, lookup, null, getConfig());
    extract.setReference(lookup);
    extract.MIN_CLIP_LENGTH = 4;
    extract.INSERT_SIZE_METRICS = new File("src/test/resources/sa.split/test.insert_size_metrics");
    extract.READ_PAIR_CONCORDANCE_METHOD = ReadPairConcordanceMethod.PERCENTAGE;
    extract.OUTPUT = output;
    extract.INPUT = input;
    try (SamReader reader = SamReaderFactory.make().open(input)) {
        extract.setup(reader.getFileHeader(), input);
    }
    List<SAMRecord> records = getRecords(input);
    List<Boolean> result = records.stream().map(r -> extract.shouldExtract(ImmutableList.of(r), lookup)[0]).collect(Collectors.toList());
    /*
		for (int i = 0; i < records.size(); i++) {
			SAMRecord r = records.get(i);
			System.out.print(r.getSupplementaryAlignmentFlag() ? "S" : " ");
			System.out.print(r.getFirstOfPairFlag() ? "1" : "2");
			System.out.print(" Extracted=" + (result.get(i) ? "y" : "n"));
			System.out.print(" HasConcPair=" + (ExtractSVReads.hasReadPairingConsistentWithReference(extract.getReadPairConcordanceCalculator(), ImmutableList.of(r)) ? "y" : "n"));
			boolean[] ra = ExtractSVReads.hasReadAlignmentConsistentWithReference(ImmutableList.of(r));
			System.out.print(" HasConcRead=" + (ra[0] ? "y" : "n") + (ra[1] ? "y" : "n"));
			System.out.println(" " + new ChimericAlignment(r).toString());
		}*/
    // primary read pair alignment implies an unexpected library fragment size = extract them all
    Assert.assertEquals(ImmutableList.of(true, true, true, true), result);
    lookup.close();
}
Also used : FixedSizeReadPairConcordanceCalculator(au.edu.wehi.idsv.FixedSizeReadPairConcordanceCalculator) Iterables(com.google.common.collect.Iterables) ReadPairConcordanceMethod(au.edu.wehi.idsv.ReadPairConcordanceMethod) Iterators(com.google.common.collect.Iterators) Lists(com.google.common.collect.Lists) ImmutableList(com.google.common.collect.ImmutableList) Files(com.google.common.io.Files) Assert.assertArrayEquals(org.junit.Assert.assertArrayEquals) StructuralVariantReadMetrics(gridss.analysis.StructuralVariantReadMetrics) ChimericAlignment(au.edu.wehi.idsv.sam.ChimericAlignment) ReferenceLookup(au.edu.wehi.idsv.picard.ReferenceLookup) IntermediateFilesTest(au.edu.wehi.idsv.IntermediateFilesTest) Assert.assertTrue(org.junit.Assert.assertTrue) IOException(java.io.IOException) Test(org.junit.Test) SynchronousReferenceLookupAdapter(au.edu.wehi.idsv.picard.SynchronousReferenceLookupAdapter) Category(org.junit.experimental.categories.Category) SamReader(htsjdk.samtools.SamReader) MetricsFile(htsjdk.samtools.metrics.MetricsFile) Collectors(java.util.stream.Collectors) File(java.io.File) SAMRecord(htsjdk.samtools.SAMRecord) List(java.util.List) Hg38Tests(au.edu.wehi.idsv.Hg38Tests) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile) Assert.assertFalse(org.junit.Assert.assertFalse) Assert(org.junit.Assert) SamReaderFactory(htsjdk.samtools.SamReaderFactory) Assert.assertEquals(org.junit.Assert.assertEquals) SamReader(htsjdk.samtools.SamReader) ReferenceLookup(au.edu.wehi.idsv.picard.ReferenceLookup) SAMRecord(htsjdk.samtools.SAMRecord) MetricsFile(htsjdk.samtools.metrics.MetricsFile) File(java.io.File) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile) IndexedFastaSequenceFile(htsjdk.samtools.reference.IndexedFastaSequenceFile) SynchronousReferenceLookupAdapter(au.edu.wehi.idsv.picard.SynchronousReferenceLookupAdapter) Category(org.junit.experimental.categories.Category)

Example 7 with ChimericAlignment

use of au.edu.wehi.idsv.sam.ChimericAlignment in project gridss by PapenfussLab.

the class SplitReadIdentificationHelper method writeSA.

/**
 * Writes the SA SAM tag for split read alignments
 * @param record primary alignment record
 * @param alignments supplementary alignments
 */
private static void writeSA(SAMRecord record, List<SAMRecord> alignments) {
    alignments.sort(ByFirstAlignedBaseReadOffset);
    List<String> satags = alignments.stream().map(r -> new ChimericAlignment(r).toString()).collect(Collectors.toList());
    record.setAttribute(SAMTag.SA.name(), satags.stream().collect(Collectors.joining(";")));
    // Primary alignment is first record as per SAMTags specs
    satags.add(0, new ChimericAlignment(record).toString());
    for (int i = 0; i < alignments.size(); i++) {
        final String currentsatag = satags.get(i + 1);
        String tag = satags.stream().filter(// can use reference equality since we're comparing to ourself
        s -> s != currentsatag).collect(Collectors.joining(";"));
        alignments.get(i).setAttribute(SAMTag.SA.name(), tag);
    }
}
Also used : Cigar(htsjdk.samtools.Cigar) Arrays(java.util.Arrays) SequenceUtil(htsjdk.samtools.util.SequenceUtil) SAMUtils(htsjdk.samtools.SAMUtils) CigarElement(htsjdk.samtools.CigarElement) CigarOperator(htsjdk.samtools.CigarOperator) ArrayUtils(org.apache.commons.lang3.ArrayUtils) SAMRecordUtil(au.edu.wehi.idsv.sam.SAMRecordUtil) Collectors(java.util.stream.Collectors) StandardCharsets(java.nio.charset.StandardCharsets) SAMRecord(htsjdk.samtools.SAMRecord) ArrayList(java.util.ArrayList) SAMTagAndValue(htsjdk.samtools.SAMRecord.SAMTagAndValue) FastqRecord(htsjdk.samtools.fastq.FastqRecord) SAMTag(htsjdk.samtools.SAMTag) List(java.util.List) Lists(com.google.common.collect.Lists) CigarUtil(au.edu.wehi.idsv.sam.CigarUtil) Comparator(java.util.Comparator) Collections(java.util.Collections) ChimericAlignment(au.edu.wehi.idsv.sam.ChimericAlignment) ChimericAlignment(au.edu.wehi.idsv.sam.ChimericAlignment)

Example 8 with ChimericAlignment

use of au.edu.wehi.idsv.sam.ChimericAlignment in project gridss by PapenfussLab.

the class SplitReadEvidence method isReference.

/**
 * Evidence provides support for no structural variant call
 * This can be due to sequence homology allowing the alignment
 * at either breakend to be placed on the opposite breakend
 */
@Override
public boolean isReference() {
    if (!isBreakendExact())
        return false;
    if (getUntemplatedSequence().length() > 0)
        return false;
    BreakendSummary location = getBreakendSummary();
    int anchorHomLen = location.nominal - location.start;
    int remoteHomLen = location.end - location.nominal;
    if (location.direction == BreakendDirection.Backward) {
        int tmp = anchorHomLen;
        anchorHomLen = remoteHomLen;
        remoteHomLen = tmp;
    }
    ChimericAlignment localAlignment = new ChimericAlignment(getSAMRecord());
    int localBases = localAlignment.getLastAlignedBaseReadOffset() - localAlignment.getFirstAlignedBaseReadOffset() + 1;
    int remoteBases = remoteAlignment.getLastAlignedBaseReadOffset() - remoteAlignment.getFirstAlignedBaseReadOffset() + 1;
    return anchorHomLen >= localBases || remoteHomLen >= remoteBases;
}
Also used : ChimericAlignment(au.edu.wehi.idsv.sam.ChimericAlignment)

Example 9 with ChimericAlignment

use of au.edu.wehi.idsv.sam.ChimericAlignment in project gridss by PapenfussLab.

the class SplitReadEvidence method create.

public static List<SplitReadEvidence> create(SAMEvidenceSource source, SAMRecord record) {
    if (record.getReadUnmappedFlag() || record.getCigar() == null)
        return Collections.emptyList();
    List<ChimericAlignment> aln = ChimericAlignment.getChimericAlignments(record);
    if (aln.isEmpty())
        return Collections.emptyList();
    if (record.getCigar().getFirstCigarElement().getOperator() == CigarOperator.HARD_CLIP || record.getCigar().getLastCigarElement().getOperator() == CigarOperator.HARD_CLIP) {
        if (!MessageThrottler.Current.shouldSupress(log, "hard clipped bases")) {
            log.warn(String.format("Read %s is hard clipped. Please run %s to soften hard clips.", record.getReadName(), ComputeSamTags.class.getName()));
        }
        record = SAMRecordUtil.hardClipToN(record);
    }
    List<SplitReadEvidence> list = new ArrayList<>(2);
    ChimericAlignment chim = new ChimericAlignment(record);
    int offset = SAMRecordUtil.getFirstAlignedBaseReadOffset(record);
    ChimericAlignment pre = aln.stream().filter(ca -> ca.getFirstAlignedBaseReadOffset() < offset).max(ChimericAlignment.ByReadOffset).orElse(null);
    ChimericAlignment post = aln.stream().filter(ca -> ca.getFirstAlignedBaseReadOffset() > offset).min(ChimericAlignment.ByReadOffset).orElse(null);
    SAMSequenceDictionary dict = source != null ? source.getContext().getDictionary() : record.getHeader().getSequenceDictionary();
    // Read is AAAXBBBYCCC
    // we are alignment B
    // pre is A
    // post is C
    // X,Y are unaligned bases
    final int rl = record.getReadLength();
    int startOffset = chim.getFirstAlignedBaseReadOffset();
    int endOffset = chim.getLastAlignedBaseReadOffset() + 1;
    if (pre != null) {
        int overlap = pre.getLastAlignedBaseReadOffset() + 1 - startOffset;
        BreakpointSummary bs = new BreakpointSummary(withOverlap(chim.predecessorBreakend(dict), overlap), withOverlap(pre.successorBreakend(dict), overlap));
        // ignore the actual alignment and just go out to the end of the read so we can assemble across multiple breakpoints
        int preStartOffset = 0;
        int preEndOffset = pre.getLastAlignedBaseReadOffset() + 1;
        if (record.getReadNegativeStrandFlag()) {
            list.add(new SplitReadEvidence(source, record, bs, rl - (INCLUDE_CLIPPED_ANCHORING_BASES ? record.getReadLength() : endOffset), rl - startOffset, rl - startOffset, rl - preEndOffset, rl - preEndOffset, rl - preStartOffset, pre, CigarUtil.widthOfImprecision(chim.cigar), CigarUtil.widthOfImprecision(pre.cigar)));
        } else {
            list.add(new SplitReadEvidence(source, record, bs, startOffset, INCLUDE_CLIPPED_ANCHORING_BASES ? record.getReadLength() : endOffset, preEndOffset, startOffset, preStartOffset, preEndOffset, pre, CigarUtil.widthOfImprecision(chim.cigar), CigarUtil.widthOfImprecision(pre.cigar)));
        }
    }
    if (post != null) {
        int overlap = endOffset - post.getFirstAlignedBaseReadOffset();
        BreakpointSummary bs = new BreakpointSummary(withOverlap(chim.successorBreakend(dict), overlap), withOverlap(post.predecessorBreakend(dict), overlap));
        int postStartOffset = post.getFirstAlignedBaseReadOffset();
        int postEndOffset = rl;
        if (record.getReadNegativeStrandFlag()) {
            list.add(new SplitReadEvidence(source, record, bs, rl - endOffset, rl - (INCLUDE_CLIPPED_ANCHORING_BASES ? 0 : startOffset), rl - postStartOffset, rl - endOffset, rl - postEndOffset, rl - postStartOffset, post, CigarUtil.widthOfImprecision(chim.cigar), CigarUtil.widthOfImprecision(post.cigar)));
        } else {
            list.add(new SplitReadEvidence(source, record, bs, INCLUDE_CLIPPED_ANCHORING_BASES ? 0 : startOffset, endOffset, endOffset, postStartOffset, postStartOffset, postEndOffset, post, CigarUtil.widthOfImprecision(chim.cigar), CigarUtil.widthOfImprecision(post.cigar)));
        }
    }
    return list;
}
Also used : ChimericAlignment(au.edu.wehi.idsv.sam.ChimericAlignment) ArrayList(java.util.ArrayList) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary)

Example 10 with ChimericAlignment

use of au.edu.wehi.idsv.sam.ChimericAlignment in project gridss by PapenfussLab.

the class ExtractSVReads method primaryAlignmentForSupplementary.

private static SAMRecord primaryAlignmentForSupplementary(SAMRecord r) {
    if (r.getSupplementaryAlignmentFlag()) {
        SAMRecord record = SAMRecordUtil.clone(r);
        List<ChimericAlignment> calist = ChimericAlignment.getChimericAlignments(r);
        if (calist.size() > 0) {
            ChimericAlignment ca = calist.get(0);
            record.setReferenceName(ca.rname);
            record.setAlignmentStart(ca.pos);
            record.setReadNegativeStrandFlag(ca.isNegativeStrand);
            record.setCigar(ca.cigar);
            record.setMappingQuality(ca.mapq);
            record.setReadBases(SAMRecord.NULL_SEQUENCE);
            record.setBaseQualities(SAMRecord.NULL_QUALS);
            return record;
        }
    }
    return r;
}
Also used : ChimericAlignment(au.edu.wehi.idsv.sam.ChimericAlignment) SAMRecord(htsjdk.samtools.SAMRecord)

Aggregations

ChimericAlignment (au.edu.wehi.idsv.sam.ChimericAlignment)20 SAMRecord (htsjdk.samtools.SAMRecord)16 Test (org.junit.Test)11 List (java.util.List)5 IntermediateFilesTest (au.edu.wehi.idsv.IntermediateFilesTest)4 SynchronousReferenceLookupAdapter (au.edu.wehi.idsv.picard.SynchronousReferenceLookupAdapter)4 ImmutableList (com.google.common.collect.ImmutableList)4 IndexedFastaSequenceFile (htsjdk.samtools.reference.IndexedFastaSequenceFile)4 File (java.io.File)4 IOException (java.io.IOException)4 ArrayList (java.util.ArrayList)4 Assert (org.junit.Assert)4 Assert.assertEquals (org.junit.Assert.assertEquals)4 Assert.assertFalse (org.junit.Assert.assertFalse)4 Assert.assertTrue (org.junit.Assert.assertTrue)4 Category (org.junit.experimental.categories.Category)4 FixedSizeReadPairConcordanceCalculator (au.edu.wehi.idsv.FixedSizeReadPairConcordanceCalculator)3 SAMRecordUtil (au.edu.wehi.idsv.sam.SAMRecordUtil)3 Lists (com.google.common.collect.Lists)3 SAMFileHeader (htsjdk.samtools.SAMFileHeader)3