SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Threshold quality score to determine the quality read of ILLUMINA reads problem edge Illumina/Solexa 35 11-02-2015 10:31 AM
Threshold quality score to determine the quality read of ILLUMINA reads problem edge General 1 09-13-2010 02:22 PM
threshold alperyilmaz RNA Sequencing 1 02-09-2010 10:09 AM
Joining Contigs coldturkey Bioinformatics 4 04-22-2009 07:58 AM
Quality score threshold? chris Bioinformatics 8 04-29-2008 12:43 AM

Reply
 
Thread Tools
Old 04-05-2012, 06:14 PM   #1
sterding
Member
 
Location: Boston

Join Date: Sep 2010
Posts: 36
Default gap threshold in joining the reads in cufflinks?

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.
Attached Images
File Type: png Screen shot 2012-04-05 at 10.17.47 PM.png (20.2 KB, 20 views)

Last edited by sterding; 04-05-2012 at 06:20 PM.
sterding is offline   Reply With Quote
Old 04-06-2012, 11:21 AM   #2
sterding
Member
 
Location: Boston

Join Date: Sep 2010
Posts: 36
Default

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.

Code:
[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;
src/common.cpp:int olap_radius = 50;
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:


Code:
[dongx@hpcc01 cufflinks-1.3.0]$ grep -P "radius" -I src/*
src/assemble.cpp:		fill_gaps(scaffolds, 2 * olap_radius);
src/assemble.cpp:		fill_gaps(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},
src/cufflinks.cpp:    fprintf(stderr, "  --overlap-radius             radius size to define overlapping reads               [ default:     50 ]\n");
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.
sterding is offline   Reply With Quote
Old 05-15-2012, 06:21 AM   #3
CHLin
Junior Member
 
Location: Taiwan

Join Date: May 2012
Posts: 3
Default 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.
CHLin is offline   Reply With Quote
Old 05-15-2012, 07:43 AM   #4
sterding
Member
 
Location: Boston

Join Date: Sep 2010
Posts: 36
Default

Quote:
Originally Posted by CHLin View Post
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.
sterding is offline   Reply With Quote
Reply

Tags
cufflinks, de novo assembly

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:45 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO