gap threshold in joining the reads in cufflinks?

[sterding]

I have a question to the de novo assembly using cufflinks: is there a threshold to join two non-overlapping reads?

As you can see from the attached figure, the top track is for alignment (BAM file from Tophat for non-strand-specific RNAseq), and the lower track is assembly from cufflinks. I run this simply as

cufflinks input.bam

It's for sure that overlapped reads are joined together into one isoform (or transfrag). What's strange for me is, why some separated reads (with gap between) are also joined, but some are not. As you can see, there are 4 gaps in the top track, where the first and last gap are not joined, while the 2nd and 3rd ones are joined. I guess it's because of the smaller gap distance in 2nd and 3rd gap, but I cannot find such 'gap threshold' in the cufflinks options.

I also tried the command with options:

cufflinks --min-intron-length 100
and
cufflinks -F 0.0 --trim-3-dropoff-frac 0 --trim-3-avgcov-thresh 0 --min-frags-per-transfrag 1 --min-intron-length 100

both don't change the result.

I appreciated if any could explain this to me. Thanks a lot.

[sterding]

OK. Without any reply yet, I dig into the source code of cufflinks by guessing something with threshold of 50. Yes, I did find some secret.

[sterding@hpcc01 cufflinks-1.3.0]$ grep -P "50\D" -I src/*
src/bundles.cpp:				if (counter.total_reads < 100 || overhang_ratio >= 0.50)
src/common.cpp:uint32_t min_intron_length = 50;
[COLOR="Red"]src/common.cpp:int olap_radius = 50;[/COLOR]
src/common.cpp:float high_phred_err_prob = 0.50; // about MAPQ = 3
src/cufflinks.cpp:    fprintf(stderr, "  --min-intron-length          minimum intron size allowed in genome                 [ default:     50 ]\n");
src/cufflinks.cpp:    fprintf(stderr, "  --overlap-radius             radius size to define overlapping reads               [ default:     50 ]\n");
src/GArgs.cpp:      <string>{;|=} e.g. disable-test;PID=S= for --disable-test PID=50 (or --PID 50) S=3.5 etc.
src/GArgs.h:       <string>{;|=} e.g. disable-test;PID=S= for --disable-test PID=50 (or --PID 50) S=3.5 etc.

See the read part: overlap_radius = 50. That's what I guessed. And this parameter seems not changeable for the user:

[dongx@hpcc01 cufflinks-1.3.0]$ grep -P "radius" -I src/*
src/assemble.cpp:		[COLOR="red"]fill_gaps[/COLOR](scaffolds, 2 * olap_radius);
src/assemble.cpp:		[COLOR="red"]fill_gaps[/COLOR](new_scaffs, 2 * olap_radius);
src/bundles.cpp:                        (!Scaffold::overlap_in_genome(L, R, olap_radius) ||
src/bundles.cpp:		if (bh->ref_id() == bundle.ref_id() && bh->left() < bundle.right() + olap_radius)
src/common.cpp:int olap_radius = 50;
src/common.h:extern int olap_radius;
src/cuffdiff.cpp:    olap_radius = 0;
src/cufflinks.cpp:{"overlap-radius",       required_argument,	     0,			 OPT_OVERLAP_RADIUS},
[COLOR="red"]src/cufflinks.cpp:    fprintf(stderr, "  --overlap-radius             radius size to define overlapping reads               [ default:     50 ]\n");[/COLOR]
src/cufflinks.cpp:				olap_radius = parseInt(0, "--overlap-radius must be at least 0", print_usage);
src/graph_optimize.cpp:        compressed.fill_gaps(2 * olap_radius);
src/scaffold_graph.cpp:				hit_j->left() - hit_i->right() < olap_radius)
src/scaffolds.cpp:								 int overlap_radius)
src/scaffolds.cpp:	int ll = lhs.left() - overlap_radius;
src/scaffolds.cpp:	int rr = rhs.right() + overlap_radius;
src/scaffolds.cpp:	int lr = lhs.right() + overlap_radius;
src/scaffolds.cpp:	int rl = rhs.left() - overlap_radius;
src/scaffolds.cpp:		if (overlap_in_genome(lhs, rhs, olap_radius))
src/scaffolds.cpp:		if (overlap_in_genome(rhs, lhs, olap_radius))
src/scaffolds.h:								  int overlap_radius);

I might re-compile the code and let the "--overlap-radius" reachable for the user. Or, It's really appreciated if cufflinks team can do this.

[CHLin]

This might help...

I think this option may help to change the parameter:
--overlap-radius [#nt]
If the gap between two units is larger than the #, they will be merged.

[sterding]

I think this option may help to change the parameter:
--overlap-radius [#nt]
If the gap between two units is larger than the #, they will be merged.

You are right; this option becomes user-settable until the recent release.