Skip to content

Commit

Permalink
modified: README.md
Browse files Browse the repository at this point in the history 
	modified:   src/TIDDIT.cpp
	modified:   src/data_structures/CoverageModule.cpp
	modified:   src/data_structures/Translocation.cpp
  • Loading branch information
@J35P312
J35P312 committed Sep 7, 2017
1 parent a467aff commit cc96413
Show file tree
Hide file tree
Showing 4 changed files with 41 additions and 82 deletions.
6 changes: 3 additions & 3 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,11 +59,11 @@ TIDDIT may be fine tuned by altering these optional parameters:
-d - the pair orientation, use this setting to override the automatic orientation selection
-p - the minimum number of discordant pairs and supplementary alignments used to call large SV. Default is 3
-p - the minimum number of discordant pairs and supplementary alignments used to call large SV. Default is 4

-r - the minimum number of supplementary alignments used to call small SV. Default is 3
-r - the minimum number of supplementary alignments used to call small SV. Default is 6
-q - the minimum mapping quality of the discordant pairs/supplementary aignments
-q - the minimum mapping quality of the discordant pairs/supplementary alignments
forming a variant. Default value is 10.
-c - the library coverage. Default is calculated from average genomic coverage.
Expand Down
6 changes: 3 additions & 3 deletions src/TIDDIT.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,8 @@ int main(int argc, char **argv) {
//MAIN VARIABLE

bool outtie = true; // library orientation
uint32_t minimumSupportingPairs = 3;
uint32_t minimumSupportingReads = 3;
uint32_t minimumSupportingPairs = 4;
uint32_t minimumSupportingReads = 6;
int min_insert = 100; // min insert size
int max_insert = 100000; // max insert size
int minimum_mapping_quality =10;
Expand All @@ -45,7 +45,7 @@ int main(int argc, char **argv) {
float insertStd;
int min_variant_size= 100;
string outputFileHeader ="output";
string version = "1.1.5";
string version = "1.1.6";

//collect all options as a vector
vector<string> arguments(argv, argv + argc);
Expand Down
23 changes: 14 additions & 9 deletions src/data_structures/CoverageModule.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,14 @@ ostream& test(ofstream &coverageOutput,string output){
//constructor
Cov::Cov(int binSize,string bamFile,string output,bool discordants){
ostream& covout=test(coverageOutput,output);
if(output != "stdout"){
coverageOutput.open((output+".tab").c_str());
static ostream& covout = coverageOutput;

}else{
static ostream& covout=cout;
}



//initialize the function
Expand All @@ -37,15 +45,8 @@ Cov::Cov(int binSize,string bamFile,string output,bool discordants){
this -> contigsNumber = 0;
this -> bamFile = bamFile;
this -> discordants=discordants;
this -> output = output;

if(output != "stdout"){
coverageOutput.open((output+".tab").c_str());
static ostream& covout = coverageOutput;

}else{
static ostream& covout=cout;
}

if (this -> discordants == true){
covout << "#CHR" << "\t" << "start" << "\t" << "end" << "\t" << "coverage" <<"\t" << "quality" << "\t" << "discordants" << endl;
}else{
Expand Down Expand Up @@ -136,7 +137,11 @@ void Cov::bin(BamAlignment currentRead){

//prints the results
void Cov::printCoverage(){
ostream& covout = test(coverageOutput,output);
ostream& covout=test(coverageOutput,this -> output);
if ( this -> output == "stdout") {
static ostream& covout=cout;
}

for(int i=0;i<contigsNumber;i++){
for(int j=0;j<coverageStructure[i].size();j++){
int binStart = j*binSize;
Expand Down
88 changes: 21 additions & 67 deletions src/data_structures/Translocation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,23 +21,13 @@ vector<string> Window::classification(int chr, int startA,int endA,double covA,i
string end= int2str(startB);

string GT="./1";
int CN=1;

double coverage= this -> meanCoverage;
float coverageToleranceFraction = 0.5/double(this -> ploidy);
float coverageTolerance=this -> meanCoverage*coverageToleranceFraction;


// for heterozygous events
CN = round(covAB / (coverage/double(this->ploidy)) ) - 1; // CN sounds like ratio to coverage per chr, minus 1? at least for het vars with the other allele at count 1.
if (CN < 0) { CN = 0; }
else if ( (CN+1 < 1.5 + 0.5*coverageToleranceFraction) && (CN+1 > 1 + 0.5*coverageToleranceFraction)) {
// dup genotype is 0/1
GT="0/1";
}

int CN = round( double(this -> ploidy)*(covAB/coverage) );
//if the orientation of one of the reads is normal, but the other read is inverted, the variant is an inversion
//both reads are pointing to the left
if(this -> pairOrientation == 0 or this -> pairOrientation == 3){
svType= "INV";
}
Expand All @@ -61,38 +51,25 @@ vector<string> Window::classification(int chr, int startA,int endA,double covA,i
svType= "DEL";
if (covAB < coverageTolerance) {
GT = "1/1";
// CN = 0;
} else {
GT = "0/1";
// CN = 1;
}
}//if the coverage of all regions is normal and region A and B are disjunct, the event is an insertion(mobile element)
else if( ( covAB < coverage+coverageTolerance and covAB > coverage-coverageTolerance ) and
( covB < coverage+coverageTolerance and covB > coverage-coverageTolerance ) and
( covA < coverage+coverageTolerance and covA > coverage-coverageTolerance )){
svType = "INS";
}
//if A and B are disjunct and only one of them have coverage above average, the event is an interspersed duplication
else if( (covA > coverage+coverageTolerance and covB < coverage+coverageTolerance ) or (covB > coverage+coverageTolerance and covA < coverage+coverageTolerance ) ){
svType = "IDUP";
// it would be a good idea to split this into one DUP and one BND, indicating the insertion point
//label the region marked as high coverage as the duplications
if(covA > coverage+coverageTolerance){
start=int2str(startA);
end = int2str(endA);
}else if(covB > coverage+coverageTolerance){
start=int2str(startB);
end=int2str(endB);
}
}
}
}

//if the event is a duplication, but the exact type cannot be specified
if(svType == "BND"){
if(covA > coverage+coverageTolerance or covB > coverage+coverageTolerance or covAB > coverage+coverageTolerance){
if(covAB > coverage+coverageTolerance){
svType = "DUP";
}
}
if (svType == "DUP" or svType == "TDUP"){
if(covAB > (double(this -> ploidy)+1)*coverage/double(this -> ploidy)+coverageTolerance){
GT = "1/1";
}else{
GT = "0/1";
}
}
vector<string> svVector;
svVector.push_back(svType);svVector.push_back(start);svVector.push_back(end);svVector.push_back(GT);svVector.push_back(int2str(CN));
return(svVector);
Expand Down Expand Up @@ -172,8 +149,8 @@ string Window::VCFHeader(string libraryData){
//set format
headerString+="##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n";
headerString+="##FORMAT=<ID=CN,Number=1,Type=Integer,Description=\"Copy number genotype for imprecise events\">\n";
headerString+="##FORMAT=<ID=PE,Number=1,Type=Integer,Description=\"Number of paired-ends that support the event\">\n";
headerString+="##FORMAT=<ID=SR,Number=1,Type=Integer,Description=\"Number of split reads that support the event\">\n";
headerString+="##FORMAT=<ID=DV,Number=1,Type=Integer,Description=\"Number of paired-ends that support the event\">\n";
headerString+="##FORMAT=<ID=RV,Number=1,Type=Integer,Description=\"Number of split reads that support the event\">\n";
//this sting contains info such as insert size and mean coverage
headerString+=libraryData+"\n";
//Header
Expand Down Expand Up @@ -1028,7 +1005,11 @@ void Window::VCFLine(map<string,float> statistics_floats,map<string,int> discord
var << this -> position2contig[chrA] << "\t" << posA << "\tSV_" << this -> numberOfEvents << "_1" << "\t" ;
var << "N" << "\t" << "<" << svType << ">";
var << "\t.\t" << filter << "\t" << infoField;
var << "\tGT:CN:PE:SR\t" << GT << ":" << CN << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";
if (CN != int2str(this -> ploidy) ){
var << "\tGT:CN:DV:RV\t" << GT << ":" << CN << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";
}else{
var << "\tGT:DV:RV\t" << GT << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";
}

SV_calls[chrA].push_back(var.str());
vector<int> row;
Expand All @@ -1041,7 +1022,7 @@ void Window::VCFLine(map<string,float> statistics_floats,map<string,int> discord
bnd_st << position2contig[chrA] << "\t" << posA << "\tSV_" << this -> numberOfEvents << "_1" << "\t";
bnd_st << "N" << "\t" << "N[" << position2contig[chrB] << ":" << posB << "[";
bnd_st << "\t.\t" << filter << "\t" << infoField;
bnd_st << "\tGT:CN:PE:SR\t" << GT << ":" << CN << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";
bnd_st << "\tGT:DV:RV\t" << GT << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";

SV_calls[chrA].push_back(bnd_st.str());
vector<int> row_st;
Expand All @@ -1055,7 +1036,7 @@ void Window::VCFLine(map<string,float> statistics_floats,map<string,int> discord
bnd_nd << position2contig[chrB] << "\t" << posB << "\tSV_" << this -> numberOfEvents << "_2" << "\t";
bnd_nd << "N" << "\t" << "N]" << position2contig[chrA] << ":" << posA << "]";
bnd_nd << "\t.\t" << filter << "\t" << infoField;
bnd_nd << "\tGT:CN:PE:SR\t" << GT << ":" << CN << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";
bnd_nd << "\tGT:DV:RV\t" << GT << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";

SV_calls[chrB].push_back(bnd_nd.str());
vector<int> row_nd;
Expand All @@ -1064,41 +1045,14 @@ void Window::VCFLine(map<string,float> statistics_floats,map<string,int> discord
SV_positions[chrB].push_back(row_nd);

}
if(svType == "IDUP"){
std::stringstream bnd_st;
bnd_st << position2contig[chrA] << "\t" << posA << "\tSV_" << this -> numberOfEvents << "_2" << "\t";
bnd_st << "N" << "\t" << "N[" << position2contig[chrB] << ":" << posB << "[";
bnd_st << "\t.\t" << filter << "\t" << infoField;
bnd_st << "\tGT:CN:PE:SR\t" << GT << ":" << CN << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";

SV_calls[chrA].push_back(bnd_st.str());
vector<int> row_st;
row_st.push_back(SV_calls[chrA].size()-1);
row_st.push_back(posA);
SV_positions[chrA].push_back(row_st);

//print the second breakend
std::stringstream bnd_nd;
bnd_nd << position2contig[chrB]<< "\t" << posB << "\tSV_" << this -> numberOfEvents << "_3" << "\t";
bnd_nd << "N" << "\t" << "N]" << position2contig[chrA] << ":" << posA << "]";
bnd_nd << "\t.\t" << filter << "\t" << infoField;
bnd_nd << "\tGT:CN:PE:SR\t" << GT << ":" << CN << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";

SV_calls[chrB].push_back(bnd_nd.str());
vector<int> row_nd;
row_nd.push_back(SV_calls[chrB].size()-1);
row_nd.push_back(posB);
SV_positions[chrB].push_back(row_nd);
}


} else {
//print the first breakend
std::stringstream bnd_st;
bnd_st << position2contig[chrA] << "\t" << posA << "\tSV_" << this -> numberOfEvents << "_1" << "\t";
bnd_st << "N" << "\t" << "N[" << position2contig[chrB] << ":" << posB << "[";
bnd_st << "\t.\t" << filter << "\t" << infoField;
bnd_st << "\tGT:CN:PE:SR\t" << GT << ":" << CN << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";
bnd_st << "\tGT:DV:RV\t" << GT << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";

SV_calls[chrA].push_back(bnd_st.str());
vector<int> row_st;
Expand All @@ -1111,7 +1065,7 @@ void Window::VCFLine(map<string,float> statistics_floats,map<string,int> discord
bnd_nd << position2contig[chrB]<< "\t" << posB << "\tSV_" << this -> numberOfEvents << "_2" << "\t";
bnd_nd << "N" << "\t" << "N]" << position2contig[chrA] << ":" << posA << "]";
bnd_nd << "\t.\t" << filter << "\t" << infoField;
bnd_nd << "\tGT:CN:PE:SR\t" << GT << ":" << CN << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";
bnd_nd << "\tGT:DV:RV\t" << GT << ":" << discordantPairStatistics["links_event"] << ":" << splitReadStatistics["split_reads"] << "\n";

SV_calls[chrB].push_back(bnd_nd.str());
vector<int> row_nd;
Expand Down

0 comments on commit cc96413

Please sign in to comment.