Skip to content

Commit

Permalink
More complete VCF output: calls that do not pass CLOVE's quality control
Browse files Browse the repository at this point in the history 
are sstill included, but with a "FAIL" in the filter field.
  • Loading branch information
@jibsch
jibsch committed Aug 28, 2015
1 parent 68bba72 commit 9cb5f62
Show file tree
Hide file tree
Showing 4 changed files with 75 additions and 32 deletions.
26 changes: 26 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,29 @@ Clove is a program designed to post-process the output of sequencing data struct
Clove scans the output of one or more sets of fusion calls for patterns of complex variations (such as balanced translocations) and summarizes fusion calls as such.
Clove also annotates the variants with read depth information from the alignment file.

----------------
INSTALLING CLOVE
----------------
Copy the latest release from the github page.
The ...-jar-with-depdendcies.jar should run without any additional requirements

-------------
RUNNING CLOVE
-------------
Invoke CLOVE with java -jar <release.jar>
This prompts you with the follwoing parameters:
Options (all mandatory -- input can be specified more than once):
-i <list of breakpoints> <algorithm (Socrates/Delly/Crest/Gustaf/BEDPE)>
-b <BAM file>
-c <mean coverage> <coverage>
-o <output filename> [default: CLOVE.vcf]
An example run of CLOVE could look like this:
java -jar clove-0.11-jar-with-dependencies.jar -i my_results.txt socrates -b my_bam.bam -c 30 7 -o my_calls.vcf
This will take the input in my_results.txt and my_bam.bam to produce calls in my_calls.vcf.

-----------
A FEW NOTES
-----------
1. The input bam file has to be sorted and indexed, as CLOVE is random accessing it.
2. The output VCF file distinguishes calls that have been classified correctly and/or pass the read depth check as "PASS" (or whatever has been provided by the original SV caller) in the filter field. All calls that failed these criteria are indicated with the "FAIL" filter. If you are interested in the "CLOVE approved" calls only, filter for anything that has not "FAIL" in the VCF entry.
3. When constructing the graph of coordinates and fusions, CLOVE discards redundant events (fusions that connect the same two nodes with identical SV type). Therefore, the output vcf is not necessarily complete with respect to the set of inputs. The algorithm will report "Events merged: X" on the command line to indicate if this has happened (X>0).
4 changes: 2 additions & 2 deletions pom.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
<groupId>au.edu.wehi</groupId>
<artifactId>clove</artifactId>
<packaging>jar</packaging>
<version>0.1</version>
<version>0.11</version>
<name>clove</name>
<url>http://github.com/PapenfussLab/clove</url>
<properties>
Expand All @@ -26,7 +26,7 @@
<configuration>
<archive>
<manifest>
<mainClass>fully.qualified.MainClass</mainClass>
<mainClass>au.edu.wehi.clove.Clove</mainClass>
</manifest>
</archive>
<descriptorRefs>
Expand Down
67 changes: 40 additions & 27 deletions src/main/java/au/edu/wehi/clove/Clove.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -505,9 +505,10 @@ public static void createVCFHeader(PrintWriter output) throws IOException{

//ALT
output.write("##ALT=<ID=DEL,Description=\"Deletion\">\n");
output.write("##ALT=<ID=DUP,Description=\"Duplication\">\n");
output.write("##ALT=<ID=TAN,Description=\"Tandem Duplication\">\n");
output.write("##ALT=<ID=INV,Description=\"Inversion\">\n");
output.write("##ALT=<ID=INS,Description=\"Insertion\">\n");
output.write("##ALT=<ID=DUP,Description=\"Complex Duplication\">\n");
output.write("##ALT=<ID=TRA,Description=\"Complex Translocation\">\n");
output.write("##ALT=<ID=CIV,Description=\"Complex Inversion\">\n");
output.write("##ALT=<ID=CVT,Description=\"Complex Inverted Translocation\">\n");
Expand Down Expand Up @@ -610,9 +611,9 @@ public static void main(String[] args) throws IOException {
continue;
Event e;
switch(algorithm){
case SOCRATES: e = Event.createNewEventFromSocratesOutput(line); break;
case SOCRATES: e = Event.createNewEventFromSocratesOutputLatest(line, count++); break;
case DELLY: e = Event.createNewEventFromDellyOutputLatest(line);break;
case CREST: e = Event.createNewEventFromCrestOutput(line); break;
case CREST: e = Event.createNewEventFromCrestOutputLatest(line, count++); break;
case GUSTAF: e = Event.createNewEventFromGustafOutput(line); if(e.size()<50) continue; break;
case BEDPE: e = Event.createNewEventFromBEDPE(line); break;
default: e = null;
Expand Down Expand Up @@ -1184,11 +1185,11 @@ else if(e2.getType() == EVENT_TYPE.INV2 ){
//this time for output
int totalEvents = 0;
for(Entry<String, TreeSet<GenomicNode>> tableEntry: genomicNodes.entrySet()) {
System.out.println("Working on Entry: "+tableEntry.toString());
//System.out.println("Working on Entry: "+tableEntry.toString());
for(GenomicNode currentNode: tableEntry.getValue()){
if(currentNode.getEvents().size() > 1){
System.out.println("Node might be shifty: "+currentNode.getEvents().size()+" members!");
System.out.println(currentNode.getEvents().get(0)+" "+currentNode.getEvents().get(1));
// System.out.println("Node might be shifty: "+currentNode.getEvents().size()+" members!");
// System.out.println(currentNode.getEvents().get(0)+" "+currentNode.getEvents().get(1));
}
totalEvents += currentNode.getEvents().size();
HashSet<Event> skipEvents = new HashSet<Event>(), deleteEvents = new HashSet<Event>(), newEvents = new HashSet<Event>();
Expand All @@ -1201,7 +1202,7 @@ else if(e2.getType() == EVENT_TYPE.INV2 ){
case INV1:
case INV2:
skipEvents.add(e);
deleteEvents.add(e);
//deleteEvents.add(e);
if(classifySimpleInversion) {
ComplexEvent e2 = new ComplexEvent(e.getC1(), e.getC2(), EVENT_TYPE.COMPLEX_INVERSION, new Event[] {e}, currentNode);
e = e2;
Expand All @@ -1213,42 +1214,39 @@ else if(e2.getType() == EVENT_TYPE.INV2 ){
e2.setQual(e2.getQual());
e2.setInfo(e2.getInfo());
}
else
continue;
else {
e.setAlt("<INV>");
e.setFailFilter();
}
break;
case DEL:
//check for deletion
//double readDepth = meanReadDepth(reader, e.getC1().getPos()+1, e.getC2().getPos()-1);
double readDepth = getReadDepth(samReader, e.getC1().getChr(), e.getC1().getPos()+1, e.getC2().getPos()-1);
skipEvents.add(e);
if(readDepth < 0 || readDepth > mean-interval){
deleteEvents.add(e);
skipEvents.add(e);
e.setId(e.getId());
e.setAlt("<DEL>");
e.setFilter(e.getFilter());
e.setQual(e.getQual());
e.setInfo(e.getInfo()+"; DP="+readDepth );
continue;
//deleteEvents.add(e);
e.setFailFilter();
} else {
System.out.print("read depth for event: "+readDepth+"\t");
//System.out.print("read depth for event: "+readDepth+"\t");
}
e.setAlt("<DEL>");
e.setInfo(e.getInfo()+"; ADP="+readDepth );
break;
case TAN:
//double readDepth = meanReadDepth(reader, e.getC1().getPos()+1, e.getC2().getPos()-1);
readDepth = getReadDepth(samReader, e.getC1().getChr(), e.getC1().getPos(), e.getC2().getPos());
skipEvents.add(e);
// //double flank = (meanReadDepth(reader, e.getC1().getPos()-200, e.getC1().getPos()) + meanReadDepth(reader, e.getC2().getPos(), e.getC2().getPos()+200))/2;
if( readDepth < 0 || readDepth < mean+interval){
//System.out.println("\t\t\t\t\t\tNot proper duplication!!");
deleteEvents.add(e);
e.setId(e.getId());
e.setAlt("<DUP>");
e.setFilter(e.getFilter());
e.setQual(e.getQual());
e.setInfo(e.getInfo()+"; DP="+readDepth );
continue;
//deleteEvents.add(e);
e.setFailFilter();
} else {
System.out.print("read depth for event: "+readDepth+"\t");
//System.out.print("read depth for event: "+readDepth+"\t");
}
e.setAlt("<DUP>");
e.setInfo(e.getInfo()+"; ADP="+readDepth );
break;
case COMPLEX_DUPLICATION:
case COMPLEX_INVERTED_DUPLICATION:
Expand All @@ -1267,9 +1265,17 @@ else if(e2.getType() == EVENT_TYPE.INV2 ){
// System.out.print("read depth for event: "+readDepth+"\t");
// }
break;
case ITX1:
case ITX2:
case INVTX1:
case INVTX2:
e.setFailFilter();
e.setAlt(Event.altVCF(e.getType()));
skipEvents.add(e);
break;
}

System.out.println(e);
//System.out.println(e);
//}
if(e.otherNode(currentNode) == currentNode){
skipEvents.add(e);
Expand All @@ -1296,12 +1302,19 @@ else if(e2.getType() == EVENT_TYPE.INV2 ){
count=0;
//reportEventComposition(genomicNodes);
/*VCF Output*/
HashSet<Event> skipEvents = new HashSet<Event>();
for(Entry<String, TreeSet<GenomicNode>> tableEntry: genomicNodes.entrySet()) {
for(GenomicNode currentNode: tableEntry.getValue()){
for(Event e: currentNode.getEvents()){
if(skipEvents.contains(e)){
continue;
} else {
skipEvents.add(e);
}
try{
writer.write(e.toVcf()+"\n");
} catch (NullPointerException npe) {
System.out.println("VCF fail");
}
}
}
Expand Down
10 changes: 7 additions & 3 deletions src/main/java/au/edu/wehi/clove/Event.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ public static Event createNewEventFromSocratesOutputLatest(String output, int co
String id="SOC"+Integer.toString(count);
String ref=".";
String qual=".";
String filter="";
String filter="PASS";

GenomicCoordinate c1 = new GenomicCoordinate(chr1, p1);
GenomicCoordinate c2 = new GenomicCoordinate(chr2, p2);
Expand Down Expand Up @@ -263,7 +263,7 @@ public static Event createNewEventFromBEDPE (String output){
String ref="";
String alt="";
String info=bits[10];
String filter = "";
String filter = "PASS";
GenomicCoordinate c1 = new GenomicCoordinate(chr1, p1);
GenomicCoordinate c2 = new GenomicCoordinate(chr2, p2);
EVENT_TYPE type = classifySocratesBreakpoint(c1, o1, c2, o2);
Expand Down Expand Up @@ -308,7 +308,7 @@ public static Event createNewEventFromCrestOutputLatest(String output, int count
String id="CRT"+Integer.toString(count);
String ref=".";
String qual=".";
String filter="";
String filter="PASS";

GenomicCoordinate c1 = new GenomicCoordinate(chr1, p1);
GenomicCoordinate c2 = new GenomicCoordinate(chr2, p2);
Expand Down Expand Up @@ -587,4 +587,8 @@ public static String altVCF(EVENT_TYPE type){
public String toVcf() {
return this.getCoord().getChr()+"\t"+this.getCoord().getPos()+"\t"+this.getId()+"\t"+this.getRef()+"\t"+this.getAlt()+"\t"+this.getQual()+"\t"+this.getFilter()+"\t"+this.getInfo();
}

public void setFailFilter(){
this.filter = "FAIL";
}
}

0 comments on commit 9cb5f62

Please sign in to comment.