Skip to content

Qsnp track unfiltered cs #388

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Mar 19, 2025
Merged

Qsnp track unfiltered cs #388

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0b87896
fix(qsnp): keep a tally of failed filter reads for compound snps
holmeso Mar 3, 2025
c06b272
refactor: removed unused code
holmeso Mar 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 50 additions & 0 deletions qannotate/test/au/edu/qimr/qannotate/modes/ConfidenceModeTest.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -294,6 +294,10 @@ public void checkMIUN() {
ConfidenceMode.checkMIUN(new String[]{"A"}, 90, "A3;C8", sb, 2, 3, null);
assertEquals("", sb.toString());

sb = new StringBuilder();
ConfidenceMode.checkMIUN(new String[]{"AA"}, 90, "AA13;AC8", sb, 2, 3, null);
assertEquals("MIUN", sb.toString());

//
}

Expand Down Expand Up @@ -811,6 +815,28 @@ public void realLifeMIUN2() {
assertEquals("PASS", vcf.getSampleFormatRecord(2).getField(VcfHeaderUtils.FORMAT_FILTER));
assertEquals("PASS", vcf.getSampleFormatRecord(3).getField(VcfHeaderUtils.FORMAT_FILTER));
assertEquals("PASS", vcf.getSampleFormatRecord(4).getField(VcfHeaderUtils.FORMAT_FILTER));

}

@Test
public void compoundSnpMIUN() {
/*
chr1 205931 . AA GG . . . GT:AD:DP:FF:FT:INF:NNS:OABS 0/0:6,0:6:AA47;AT1;A_1;GG3:.:.:.:AA3[]3[] 0/1:20,4:25:AA182;GG2;TA1:.:SOMATIC:4:AA
7[]13[];GA1[]0[];GG0[]4[]
*/
VcfRecord vcf = VcfUtils.createVcfRecord(ChrPositionUtils.getChrPosition("chr1", 205931, 205931), ".", "CA", "GG");
vcf.setInfo("FLANK=GTAAAACTGGA;BaseQRankSum=0.325;ClippingRankSum=0.000;DP=58;ExcessHet=3.0103;FS=4.683;MQ=55.10;MQRankSum=-6.669;QD=4.63;ReadPosRankSum=-0.352;SOR=1.425;IN=1;DB;VLD;HOM=3,TATATGTAAAgCTGGATTAAT;EFF=downstream_gene_variant(MODIFIER||914|||MST1P2|unprocessed_pseudogene|NON_CODING|ENST00000457982||1),intergenic_region(MODIFIER||||||||||1)");
vcf.setFormatFields(java.util.Arrays.asList(
"GT:AD:DP:FF:FT:INF:NNS:OABS",
"0/0:36,0:36:CA12;GG4;_A2:.:.:.:CA17[]19[];C_1[]0[]",
"0/1:102,11:114:AA1;CA30;CC1;CT1;C_3;GG50;GT1:.:SOMATIC:10:AA1[]0[];CA61[]41[];GG4[]7[];G_1[]0[];_A0[]2[]"));
ConfidenceMode cm = new ConfidenceMode(TWO_SAMPLE_ONE_CALLER_META);
cm.positionRecordMap.put(vcf.getChrPosition(), List.of(vcf));
cm.addAnnotation();
vcf = cm.positionRecordMap.get(vcf.getChrPosition()).getFirst();
assertEquals("MIUN", vcf.getSampleFormatRecord(1).getField(VcfHeaderUtils.FORMAT_FILTER));
assertEquals("PASS", vcf.getSampleFormatRecord(2).getField(VcfHeaderUtils.FORMAT_FILTER));

}

@Test
Expand Down Expand Up @@ -1345,6 +1371,30 @@ public void confidenceRealLifeMerged9() {
assertEquals("PASS", vcf.getSampleFormatRecord(4).getField(VcfHeaderUtils.FORMAT_FILTER));
}

@Test
public void realLifeCSMIUN() {
/*
chr1 11445731 rs386628485 AG GC . . IN=1;DB;HOM=0,ACAGAGAGACagAGAGTCAGAG GT:AD:DP:FF:FT:INF:NNS:OABS 0/0:18,1:20:AC1;AG6;AGC1;AT1;A_1;CC1;GC7:PASS:.:.:AG7[]11[];GC0[]1[];GG0
[]1[];_C1[]0[] 0/1:32,4:36:AA1;AG11;A_1;GA1;GC18;G_1:MR:SOMATIC:4:AG21[]11[];A_1[]0[];GC2[]2[];_C1[]0[] ./.:.:.:.:COV:.:.:. ./.:.:.:.:COV:.:.:.
*/
VcfRecord vcf = new VcfRecord(new String[]{"chr1", "11445731", "rs386628485", "AG", "GC", ".", ".", "IN=1;DB;HOM=0,ACAGAGAGACagAGAGTCAGAG"
, "GT:AD:DP:FF:FT:INF:NNS:OABS"
, "0/0:18,1:20:AC1;AG6;AGC1;AT1;A_1;CC1;GC7:.:.:.:AG7[]11[];GC0[]1[];GG0[]1[];_C1[]0[]"
, "0/1:32,4:36:AA1;AG11;A_1;GA1;GC18;G_1:.:SOMATIC:4:AG21[]11[];A_1[]0[];GC2[]2[];_C1[]0[]"
, "./.:.:.:.:.:.:.:."
, "./.:.:.:.:.:.:.:."});
ConfidenceMode cm = new ConfidenceMode(TWO_SAMPLE_TWO_CALLER_META);
cm.positionRecordMap.put(vcf.getChrPosition(), List.of(vcf));
cm.addAnnotation();

vcf = cm.positionRecordMap.get(vcf.getChrPosition()).getFirst();
assertEquals("MIUN", vcf.getSampleFormatRecord(1).getField(VcfHeaderUtils.FORMAT_FILTER));
assertEquals("MR", vcf.getSampleFormatRecord(2).getField(VcfHeaderUtils.FORMAT_FILTER));
assertEquals("COV", vcf.getSampleFormatRecord(3).getField(VcfHeaderUtils.FORMAT_FILTER));
assertEquals("COV", vcf.getSampleFormatRecord(4).getField(VcfHeaderUtils.FORMAT_FILTER));

}

@Test
public void applyMRFilter() {
assertFalse(ConfidenceMode.applyMutantReadFilter(null, null, -1));
Expand Down
57 changes: 40 additions & 17 deletions qcommon/src/org/qcmg/common/model/Accumulator.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,10 +43,10 @@ public class Accumulator {

private final int position;

private short failedFilterACount = 0;
private short failedFilterCCount = 0;
private short failedFilterGCount = 0;
private short failedFilterTCount = 0;
private TLongList failedFilterACount;
private TLongList failedFilterCCount;
private TLongList failedFilterGCount;
private TLongList failedFilterTCount;

private TLongList readNameHashStrandBasePositionQualities;

Expand All @@ -58,19 +58,23 @@ public int getPosition() {
return position;
}

public void addFailedFilterBase(final byte base) {
public void addFailedFilterBase(final byte base, long readNameHash) {
switch (base) {
case A_BYTE:
failedFilterACount++;
if (null == failedFilterACount) failedFilterACount = new TLongArrayList();
failedFilterACount.add(readNameHash);
break;
case C_BYTE:
failedFilterCCount++;
if (null == failedFilterCCount) failedFilterCCount = new TLongArrayList();
failedFilterCCount.add(readNameHash);
break;
case G_BYTE:
failedFilterGCount++;
if (null == failedFilterGCount) failedFilterGCount = new TLongArrayList();
failedFilterGCount.add(readNameHash);
break;
case T_BYTE:
failedFilterTCount++;
if (null == failedFilterTCount) failedFilterTCount = new TLongArrayList();
failedFilterTCount.add(readNameHash);
break;
default: /* do nothing */
break;
Expand Down Expand Up @@ -120,17 +124,17 @@ public String toString() {

public String getFailedFilterPileup() {
StringBuilder sb = new StringBuilder();
if (failedFilterACount > 0) {
StringUtils.updateStringBuilder(sb, A_STRING + failedFilterACount, Constants.SEMI_COLON);
if (null != failedFilterACount && ! failedFilterACount.isEmpty()) {
StringUtils.updateStringBuilder(sb, A_STRING + failedFilterACount.size(), Constants.SEMI_COLON);
}
if (failedFilterCCount > 0) {
StringUtils.updateStringBuilder(sb, C_STRING + failedFilterCCount, Constants.SEMI_COLON);
if (null != failedFilterCCount && ! failedFilterCCount.isEmpty()) {
StringUtils.updateStringBuilder(sb, C_STRING + failedFilterCCount.size(), Constants.SEMI_COLON);
}
if (failedFilterGCount > 0) {
StringUtils.updateStringBuilder(sb, G_STRING + failedFilterGCount, Constants.SEMI_COLON);
if (null != failedFilterGCount && ! failedFilterGCount.isEmpty()) {
StringUtils.updateStringBuilder(sb, G_STRING + failedFilterGCount.size(), Constants.SEMI_COLON);
}
if (failedFilterTCount > 0) {
StringUtils.updateStringBuilder(sb, T_STRING + failedFilterTCount, Constants.SEMI_COLON);
if (null != failedFilterTCount && ! failedFilterTCount.isEmpty()) {
StringUtils.updateStringBuilder(sb, T_STRING + failedFilterTCount.size(), Constants.SEMI_COLON);
}
return !sb.isEmpty() ? sb.toString() : Constants.MISSING_DATA_STRING;
}
Expand All @@ -140,4 +144,23 @@ public int getCoverage() {
return null == readNameHashStrandBasePositionQualities ? 0 : readNameHashStrandBasePositionQualities.size() / 2;
}

public boolean isEmpty() {
return null == readNameHashStrandBasePositionQualities && null == failedFilterACount && null == failedFilterCCount && null == failedFilterGCount && null == failedFilterTCount;
}

public TLongList getFailedFilterACount() {
return failedFilterACount;
}

public TLongList getFailedFilterCCount() {
return failedFilterCCount;
}

public TLongList getFailedFilterGCount() {
return failedFilterGCount;
}

public TLongList getFailedFilterTCount() {
return failedFilterTCount;
}
}
14 changes: 8 additions & 6 deletions qcommon/src/org/qcmg/common/util/AccumulatorUtils.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,6 @@ public class AccumulatorUtils {
public static final long T_BASE_BIT = 0x400000000000000L;
public static final int T_BASE_BIT_POSITION = 58;


public static final long STRAND_BIT = 0x8000000000000000L;
public static final int STRAND_BIT_POSITION = 63;
public static final long END_OF_READ_BIT = 0x4000000000000000L;
Expand All @@ -69,7 +68,7 @@ public class AccumulatorUtils {

/**
* This removes reads that have the same read name hash from the accumulator.
* If
* <p>
* If the duplicates have the same base, then 1 is left, if they have different bases, they are both (all?) removed
* <p>
* This method updates the Accumulator object that is passed in, and is therefore not side-effect free
Expand Down Expand Up @@ -236,6 +235,7 @@ public static int[] getBaseCountByStrand(int[] array, char c) {
* strand (bit 63)
* end of read (bit 62)
* base (bits 58-61)
* passedFilter (bit 57)
* quality (bits 32-40)
* position (bits 0-31)
*
Expand Down Expand Up @@ -1014,8 +1014,9 @@ public static TLongIntMap getReadNameHashStartPositionMap(Accumulator acc) {
if (null != acc) {
TLongList list = acc.getData();
if (null != list) {
TLongIntMap map = new TLongIntHashMap(list.size() * 2);
for (int i = 0, len = list.size(); i < len; i += 2) {
int len = list.size();
TLongIntMap map = new TLongIntHashMap(len);
for (int i = 0; i < len; i += 2) {

int startPosition = (int) list.get(i + 1);
if (((list.get(i + 1) >>> STRAND_BIT_POSITION) & 1) == 0) {
Expand All @@ -1038,8 +1039,9 @@ public static TLongCharMap getReadNameHashBaseMap(Accumulator acc) {
if (null != acc) {
TLongList list = acc.getData();
if (null != list) {
TLongCharMap map = new TLongCharHashMap(list.size() * 2);
for (int i = 0, len = list.size(); i < len; i += 2) {
int len = list.size();
TLongCharMap map = new TLongCharHashMap(len);
for (int i = 0; i < len; i += 2) {

char base = getBaseAsCharFromLong(list.get(i + 1));
if (((list.get(i + 1) >>> STRAND_BIT_POSITION) & 1) == 0) {
Expand Down
29 changes: 16 additions & 13 deletions qcommon/test/org/qcmg/common/model/AccumulatorTest.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,11 +81,12 @@ public void endOfReads() {
public void testUnfilteredPileup() {
Accumulator acc = new Accumulator(1);
String basesString = "ACGT";
for (byte b : basesString.getBytes()) acc.addFailedFilterBase(b);
long readNameHash = 1;
for (byte b : basesString.getBytes()) acc.addFailedFilterBase(b, readNameHash);
assertEquals("A1;C1;G1;T1", acc.getFailedFilterPileup());
for (byte b : basesString.getBytes()) acc.addFailedFilterBase(b);
for (byte b : basesString.getBytes()) acc.addFailedFilterBase(b, readNameHash + 1);
assertEquals("A2;C2;G2;T2", acc.getFailedFilterPileup());
for (byte b : basesString.getBytes()) acc.addFailedFilterBase(b);
for (byte b : basesString.getBytes()) acc.addFailedFilterBase(b, readNameHash + 2);
assertEquals("A3;C3;G3;T3", acc.getFailedFilterPileup());
}

Expand All @@ -97,16 +98,17 @@ public void testUnfilteredPileupPercentage() {
}

String basesString = "GG";
long readNameHash = 1;
for (byte b : basesString.getBytes()) {
acc.addFailedFilterBase(b);
acc.addFailedFilterBase(b, readNameHash++);
}
assertEquals("G2", acc.getFailedFilterPileup());
/*
* need 3 percent
*/
basesString = "G";
for (byte b : basesString.getBytes()) {
acc.addFailedFilterBase(b);
acc.addFailedFilterBase(b, readNameHash++);
}
assertEquals("G3", acc.getFailedFilterPileup());
}
Expand All @@ -115,31 +117,32 @@ public void testUnfilteredPileupPercentage() {
@Test
public void singleUnfilteredPileup() {
Accumulator acc = new Accumulator(1);
for (byte b : "ACGT".getBytes()) acc.addFailedFilterBase(b);
long readNameHash = 1;
for (byte b : "ACGT".getBytes()) acc.addFailedFilterBase(b, readNameHash++);
assertEquals("A1;C1;G1;T1", acc.getFailedFilterPileup());

acc = new Accumulator(1);
for (byte b : "ACGTA".getBytes()) acc.addFailedFilterBase(b);
for (byte b : "ACGTA".getBytes()) acc.addFailedFilterBase(b, readNameHash++);
assertEquals("A2;C1;G1;T1", acc.getFailedFilterPileup());

acc = new Accumulator(1);
for (byte b : "ACCGT".getBytes()) acc.addFailedFilterBase(b);
for (byte b : "ACCGT".getBytes()) acc.addFailedFilterBase(b, readNameHash++);
assertEquals("A1;C2;G1;T1", acc.getFailedFilterPileup());

acc = new Accumulator(1);
for (byte b : "ATTTGT".getBytes()) acc.addFailedFilterBase(b);
for (byte b : "ATTTGT".getBytes()) acc.addFailedFilterBase(b, readNameHash++);
assertEquals("A1;G1;T4", acc.getFailedFilterPileup());

acc = new Accumulator(1);
for (byte b : "AAAATTTGT".getBytes()) acc.addFailedFilterBase(b);
for (byte b : "AAAATTTGT".getBytes()) acc.addFailedFilterBase(b, readNameHash++);
assertEquals("A4;G1;T4", acc.getFailedFilterPileup());

acc = new Accumulator(1);
for (byte b : "AAAACTTTCGT".getBytes()) acc.addFailedFilterBase(b);
for (byte b : "AAAACTTTCGT".getBytes()) acc.addFailedFilterBase(b, readNameHash++);
assertEquals("A4;C2;G1;T4", acc.getFailedFilterPileup());

acc = new Accumulator(1);
for (byte b : "AAAACTTTCGTG".getBytes()) acc.addFailedFilterBase(b);
for (byte b : "AAAACTTTCGTG".getBytes()) acc.addFailedFilterBase(b, readNameHash++);
assertEquals("A4;C2;G2;T4", acc.getFailedFilterPileup());
}

Expand Down Expand Up @@ -217,7 +220,7 @@ public void getGenotypeRealLife() {
Accumulator acc = new Accumulator(1);
for (int i = 1; i <= 60; i++) acc.addBase((byte) 'G', (byte) 40, false, 1, 1, 2, i);
for (int i = 1; i <= 5; i++) acc.addBase((byte) 'C', (byte) 42, false, 1, 1, 2, i + 61);
for (int i = 1; i <= 1; i++) acc.addBase((byte) 'C', (byte) 42, true, 1, 1, 2, i + 67);
for (int i = 1; i == 1; i++) acc.addBase((byte) 'C', (byte) 42, true, 1, 1, 2, i + 67);

assertEquals("C1[42]5[42];G0[0]60[40]", AccumulatorUtils.getOABS(acc));
/*
Expand Down
Loading