Skip to content

Commit 1e3080c

Browse files
TimD1TimD1
committed
feat: added per-read anchor requirement to junction extract
 - closes griffithlab#186, reads now only 'support' a junction if they have at least a given minimum anchor length, supplied with the '-A' flag (default 0)
1 parent 409a492 commit 1e3080c

File tree

4 files changed

+34
-9
lines changed

4 files changed

+34
-9
lines changed

‎src/cis-splice-effects/cis_splice_effects_identifier.cc

+8-2
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,8 @@ void CisSpliceEffectsIdentifier::usage(ostream& out) {
4444
out << "\t\t" << "-t STR\tTag used in bam to label strand. [XS]" << endl;
4545
out << "\t\t" << "-a INT\tMinimum anchor length. Junctions which satisfy a minimum \n"
4646
<< "\t\t\t " << "anchor length on both sides are reported. [8]" << endl;
47+
out << "\t\t" << "-A INT\tMinimum read anchor length. Reads which satisfy a minimum \n"
48+
<< "\t\t\t " << "anchor length on both sides 'support' a junction. [0]" << endl;
4749
out << "\t\t" << "-m INT\tMinimum intron length. [70]" << endl;
4850
out << "\t\t" << "-M INT\tMaximum intron length. [500000]" << endl;
4951
out << "\t\t" << "-f INT\tOnly use alignments where all flag bits set here are set. [0]" << endl;
@@ -116,7 +118,7 @@ void CisSpliceEffectsIdentifier::parse_options(int argc, char* argv[]) {
116118
optind = 1; //Reset before parsing again.
117119
stringstream help_ss;
118120
char c;
119-
while((c = getopt(argc, argv, "o:w:v:j:e:Ei:ISht:s:a:m:M:f:F:q:b:C")) != -1) {
121+
while((c = getopt(argc, argv, "o:w:v:j:e:Ei:ISht:s:a:A:m:M:f:F:q:b:C")) != -1) {
120122
switch(c) {
121123
case 'o':
122124
output_file_ = string(optarg);
@@ -167,6 +169,9 @@ void CisSpliceEffectsIdentifier::parse_options(int argc, char* argv[]) {
167169
case 'a':
168170
min_anchor_length_ = atoi(optarg);
169171
break;
172+
case 'A':
173+
min_read_anchor_length_ = atoi(optarg);
174+
break;
170175
case 'm':
171176
min_intron_length_ = atoi(optarg);
172177
break;
@@ -298,7 +303,8 @@ void CisSpliceEffectsIdentifier::identify() {
298303
ref_to_pass = "NA";
299304
}
300305
JunctionsExtractor je1(bam_, variant_region, strandness_,
301-
strand_tag_, min_anchor_length_, min_intron_length_, max_intron_length_,
306+
strand_tag_, min_anchor_length_, min_read_anchor_length_,
307+
min_intron_length_, max_intron_length_,
302308
filter_flags_, require_flags_, min_map_qual_, ref_to_pass);
303309
je1.identify_junctions_from_BAM();
304310
vector<Junction> junctions = je1.get_all_junctions();

‎src/cis-splice-effects/cis_splice_effects_identifier.h

+4-2
Original file line numberDiff line numberDiff line change
@@ -87,9 +87,10 @@ class CisSpliceEffectsIdentifier {
8787
//tag used in BAM to denote strand, default "XS"
8888
string strand_tag_;
8989
//Minimum anchor length for junctions
90-
//Junctions need at least this many bp overlap
91-
// on both ends.
90+
//Junctions need at least this many bp overlap on both ends.
9291
uint32_t min_anchor_length_;
92+
//Reads need at least this many bp overlap on both ends to support junction.
93+
uint32_t min_read_anchor_length_;
9394
//Minimum length of an intron, i.e min junction width
9495
uint32_t min_intron_length_;
9596
//Maximum length of an intron, i.e max junction width
@@ -118,6 +119,7 @@ class CisSpliceEffectsIdentifier {
118119
strandness_(-1),
119120
strand_tag_("XS"),
120121
min_anchor_length_(8),
122+
min_read_anchor_length_(0),
121123
min_intron_length_(70),
122124
max_intron_length_(500000),
123125
filter_flags_(0),

‎src/junctions/junctions_extractor.cc

+15-2
Original file line numberDiff line numberDiff line change
@@ -43,14 +43,17 @@ int JunctionsExtractor::parse_options(int argc, char *argv[]) {
4343
optind = 1; //Reset before parsing again.
4444
int c;
4545
stringstream help_ss;
46-
while((c = getopt(argc, argv, "ha:m:M:f:F:q:o:r:t:s:b:")) != -1) {
46+
while((c = getopt(argc, argv, "ha:A:m:M:f:F:q:o:r:t:s:b:")) != -1) {
4747
switch(c) {
4848
case 'h':
4949
usage(help_ss);
5050
throw common::cmdline_help_exception(help_ss.str());
5151
case 'a':
5252
min_anchor_length_ = atoi(optarg);
5353
break;
54+
case 'A':
55+
min_read_anchor_length_ = atoi(optarg);
56+
break;
5457
case 'm':
5558
min_intron_length_ = atoi(optarg);
5659
break;
@@ -123,6 +126,7 @@ int JunctionsExtractor::parse_options(int argc, char *argv[]) {
123126
}
124127

125128
cerr << "Minimum junction anchor length: " << min_anchor_length_ << endl;
129+
cerr << "Minimum read anchor length: " << min_read_anchor_length_ << endl;
126130
cerr << "Minimum intron length: " << min_intron_length_ << endl;
127131
cerr << "Maximum intron length: " << max_intron_length_ << endl;
128132
cerr << "Require alignment flags: " << require_flags_ << endl;
@@ -144,6 +148,8 @@ int JunctionsExtractor::usage(ostream& out) {
144148
out << "Options:" << endl;
145149
out << "\t\t" << "-a INT\tMinimum anchor length. Junctions which satisfy a minimum \n"
146150
<< "\t\t\t " << "anchor length on both sides are reported. [8]" << endl;
151+
out << "\t\t" << "-A INT\tMinimum read anchor length. Reads which satisfy a minimum \n"
152+
<< "\t\t\t " << "anchor length on both sides 'support' a junction. [0]" << endl;
147153
out << "\t\t" << "-m INT\tMinimum intron length. [70]" << endl;
148154
out << "\t\t" << "-M INT\tMaximum intron length. [500000]" << endl;
149155
out << "\t\t" << "-f INT\tOnly use alignments where all flag bits set here are set. [0]" << endl;
@@ -175,12 +181,19 @@ string JunctionsExtractor::get_new_junction_name() {
175181
return name_ss.str();
176182
}
177183

178-
//Do some basic qc on the junction
184+
//Update if junction passes QC based on current read alignment
179185
bool JunctionsExtractor::junction_qc(Junction &j1) {
186+
// don't add support for junction if intron is wrong size
180187
if(j1.end - j1.start < min_intron_length_ ||
181188
j1.end - j1.start > max_intron_length_) {
182189
return false;
183190
}
191+
192+
// don't add support for junction if read isn't sufficiently anchored
193+
if(j1.start - j1.thick_start < min_read_anchor_length_) return false;
194+
if(j1.thick_end - j1.end < min_read_anchor_length_) return false;
195+
196+
// add support, update if this junction is sufficiently anchored
184197
if(j1.start - j1.thick_start >= min_anchor_length_)
185198
j1.has_left_min_anchor = true;
186199
if(j1.thick_end - j1.end >= min_anchor_length_)

‎src/junctions/junctions_extractor.h

+7-3
Original file line numberDiff line numberDiff line change
@@ -153,9 +153,10 @@ class JunctionsExtractor {
153153
//Reference FASTA file
154154
string ref_;
155155
//Minimum anchor length for junctions
156-
//Junctions need atleast this many bp overlap
157-
// on both ends.
156+
//Junctions need at least this many bp overlap on both ends.
158157
uint32_t min_anchor_length_;
158+
//Reads need at least this many bp overlap to support a junction
159+
uint32_t min_read_anchor_length_;
159160
//Minimum length of an intron, i.e min junction width
160161
uint32_t min_intron_length_;
161162
//Maximum length of an intron, i.e max junction width
@@ -190,6 +191,7 @@ class JunctionsExtractor {
190191
//Default constructor
191192
JunctionsExtractor() {
192193
min_anchor_length_ = 8;
194+
min_read_anchor_length_ = 0;
193195
min_intron_length_ = 70;
194196
max_intron_length_ = 500000;
195197
filter_flags_ = 0;
@@ -211,6 +213,7 @@ class JunctionsExtractor {
211213
int strandness1,
212214
string strand_tag1,
213215
uint32_t min_anchor_length1,
216+
uint32_t min_read_anchor_length1,
214217
uint32_t min_intron_length1,
215218
uint32_t max_intron_length1,
216219
uint16_t filter_flags,
@@ -222,7 +225,8 @@ class JunctionsExtractor {
222225
strandness_(strandness1),
223226
strand_tag_(strand_tag1),
224227
min_anchor_length_(min_anchor_length1),
225-
min_intron_length_(min_anchor_length1),
228+
min_read_anchor_length_(min_read_anchor_length1),
229+
min_intron_length_(min_intron_length1),
226230
max_intron_length_(max_intron_length1),
227231
filter_flags_(filter_flags),
228232
require_flags_(require_flags),

0 commit comments

Comments
 (0)