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();
}
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);
}
}
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;
}
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;
}
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;
}
Aggregations