No announcement yet.
  • Filter
  • Time
  • Show
Clear All
new posts

  • 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
    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 Files
    Last edited by sterding; 04-05-2012, 06:20 PM.

  • #2
    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.


    • #3
      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.


      • #4
        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.


        Latest Articles


        • seqadmin
          Advanced Methods for the Detection of Infectious Disease
          by seqadmin

          The recent pandemic caused worldwide health, economic, and social disruptions with its reverberations still felt today. A key takeaway from this event is the need for accurate and accessible tools for detecting and tracking infectious diseases. Timely identification is essential for early intervention, managing outbreaks, and preventing their spread. This article reviews several valuable tools employed in the detection and surveillance of infectious diseases.
          11-27-2023, 01:15 PM
        • seqadmin
          Strategies for Investigating the Microbiome
          by seqadmin

          Microbiome research has led to the discovery of important connections to human and environmental health. Sequencing has become a core investigational tool in microbiome research, a subject that we covered during a recent webinar. Our expert speakers shared a number of advancements including improved experimental workflows, research involving transmission dynamics, and invaluable analysis resources. This article recaps their informative presentations, offering insights...
          11-09-2023, 07:02 AM





        Topics Statistics Last Post
        Started by seqadmin, Yesterday, 08:23 AM
        0 responses
        Last Post seqadmin  
        Started by seqadmin, 12-01-2023, 09:55 AM
        0 responses
        Last Post seqadmin  
        Started by seqadmin, 11-30-2023, 10:48 AM
        0 responses
        Last Post seqadmin  
        Started by seqadmin, 11-29-2023, 08:26 AM
        0 responses
        Last Post seqadmin