Merge in Rory Kirshner's changes to cufflinks and cuffdiff. These changes add normalization support to cufflinks and cuffdiff and replicate support for cuffdiff.

2026-10-25 00:12:08 +08:00 · 2011-01-21 13:30:13 -05:00
parent 97fcb34621 ec2c3cd963
commit 6ae2ba0f09
4 changed files with 124 additions and 11 deletions
@@ -2,6 +2,33 @@

 import optparse, os, shutil, subprocess, sys, tempfile

+def group_callback( option, op_str, value, parser ):
+    groups = []
+    flist = []
+    for arg in parser.rargs:
+        arg = arg.strip()
+        if arg[0] is "-":
+            break
+        elif arg[0] is ",":
+            groups.append(flist)
+            flist = []
+        else:
+            flist.append(arg)
+    groups.append(flist)
+
+    setattr(parser.values, option.dest, groups)
+    
+def label_callback( option, op_str, value, parser ):
+    labels = []
+    for arg in parser.rargs:
+        arg = arg.strip()
+        if arg[0] is "-":
+            break
+        else:
+            labels.append(arg)
+
+    setattr(parser.values, option.dest, labels)
+
 def stop_err( msg ):
    sys.stderr.write( "%s\n" % msg )
    sys.exit()
@@ -39,10 +66,17 @@ def __main__():
    parser.add_option( '--max-mle-iterations', dest='max_mle_iterations', help='Sets the number of iterations allowed during maximum likelihood estimation of abundances. Default: 5000' )
    
    # Wrapper / Galaxy options.
+    parser.add_option( '-f', '--files', dest='groups', action="callback", callback=group_callback, help="Groups to be processed, groups are separated by spaces, replicates in a group comma separated. group1_rep1,group1_rep2 group2_rep1,group2_rep2, ..., groupN_rep1, groupN_rep2" )
    parser.add_option( '-A', '--inputA', dest='inputA', help='A transcript GTF file produced by cufflinks, cuffcompare, or other source.')
    parser.add_option( '-1', '--input1', dest='input1', help='File of RNA-Seq read alignments in the SAM format. SAM is a standard short read alignment, that allows aligners to attach custom tags to individual alignments, and Cufflinks requires that the alignments you supply have some of these tags. Please see Input formats for more details.' )
    parser.add_option( '-2', '--input2', dest='input2', help='File of RNA-Seq read alignments in the SAM format. SAM is a standard short read alignment, that allows aligners to attach custom tags to individual alignments, and Cufflinks requires that the alignments you supply have some of these tags. Please see Input formats for more details.' )
+
+    # Label options
+    parser.add_option('-L', '--labels', dest='labels', action="callback", callback=label_callback, help="Labels for the groups the replicates are in.")
    
+	# Normalization options.
+    parser.add_option( "-N", "--quartile-normalization", dest="do_normalization", action="store_true" )
+
    # Bias correction options.
    parser.add_option( '-r', dest='do_bias_correction', action="store_true", help='Providing Cufflinks with a multifasta file via this option instructs it to run our new bias detection and correction algorithm which can significantly improve accuracy of transcript abundance estimates.')
    parser.add_option( '', '--dbkey', dest='dbkey', help='The build of the reference dataset' )
@@ -121,14 +155,30 @@ def __main__():
        cmd += ( " --num-importance-samples %i" % int ( options.num_importance_samples ) )
    if options.max_mle_iterations:
        cmd += ( " --max-mle-iterations %i" % int ( options.max_mle_iterations ) )
+    if options.do_normalization:
+        cmd += ( " -N" )
    if options.do_bias_correction:
        cmd += ( " -r %s" % seq_path )
+            
+    # Add inputs.
+    # For replicate analysis: group1_rep1,group1_rep2 groupN_rep1,groupN_rep2
+    if options.groups:
+        cmd += " --labels "
+        for label in options.labels:
+            cmd += label + ","
+        cmd = cmd[:-1]
+
+        cmd += " " + options.inputA + " "
+
+        for group in options.groups:
+            for filename in group:
+                cmd += filename + ","
+            cmd = cmd[:-1] + " "
+    else: 
+        cmd += " " + options.inputA + " " + options.input1 + " " + options.input2
        
    # Debugging.
    print cmd
-            
-    # Add inputs.
-    cmd += " " + options.inputA + " " + options.input1 + " " + options.input2

    # Run command.
    try:
@@ -22,15 +22,17 @@
            --cds_diff_output=$cds_diff
            --promoters_diff_output=$promoters_diff
            
-            --inputA=$gtf_input
-            --input1=$aligned_reads1
-            --input2=$aligned_reads2
-            
            ## Set paired-end data parameters?
            #if $singlePaired.sPaired == "Yes":
                -m $singlePaired.mean_inner_distance
                -s $singlePaired.inner_distance_std_dev
+                #end if
+
+            ## Normalization?
+            #if str($do_normalization) == "Yes":
+            -N
            #end if
+
            
            ## Bias correction?
            #if $bias_correction.do_bias_correction == "Yes":
@@ -43,14 +45,56 @@
                --dbkey=${gtf_input.metadata.dbkey} 
                --index_dir=${GALAXY_DATA_INDEX_DIR}
            #end if
+                
+            ## Inputs.
+            --inputA=$gtf_input
+            #if $group_analysis.do_groups == "No":
+                --input1=$aligned_reads1
+                --input2=$aligned_reads2
+            #else:
+                ## Replicates.
+                --labels
+                #for $group in $group_analysis.groups
+                    ${group.group}
+                #end for
+                --files
+                #for $group in $group_analysis.groups
+                    #for $file in $group.files:
+                        ${file.file}
+                    #end for
+                    ,
+                #end for
+            #end if
+
    </command>
    <inputs>
        <param format="gtf" name="gtf_input" type="data" label="Transcripts" help="A transcript GTF file produced by cufflinks, cuffcompare, or other source."/>
-        <param format="sam,bam" name="aligned_reads1" type="data" label="SAM or BAM file of aligned RNA-Seq reads" help=""/>
-        <param format="sam,bam" name="aligned_reads2" type="data" label="SAM or BAM file of aligned RNA-Seq reads" help=""/>
+        <conditional name="group_analysis"> 
+            <param name="do_groups" type="select" label="Perform replicate analysis" help="Perform cuffdiff with replicates in each group.">
+                <option value="No">No</option>
+                <option value="Yes">Yes</option>
+            </param>
+            <when value="Yes">
+                <repeat name="groups" title="Group">
+                    <param name="group" title="Group name" type="text" label="Group name (no spaces or commas)"/>
+                    <repeat name="files" title="Replicate">
+                        <param name="file" label="Add file" type="data" format="sam,bam"/>
+                    </repeat>
+                </repeat>
+            </when>
+            <when value="No">
+                <param format="sam,bam" name="aligned_reads1" type="data" label="SAM or BAM file of aligned RNA-Seq reads" help=""/>
+                <param format="sam,bam" name="aligned_reads2" type="data" label="SAM or BAM file of aligned RNA-Seq reads" help=""/>
+            </when>
+        </conditional>
+
        <param name="fdr" type="float" value="0.05" label="False Discovery Rate" help="The allowed false discovery rate."/>
        <param name="min_mapqual" type="integer" value="0" label="Min SAM Mapping Quality" help="Instructs Cufflinks to ignore alignments with a SAM mapping quality lower than this number."/>
        <param name="min_alignment_count" type="integer" value="0" label="Min Alignment Count" help="The minimum number of alignments in a locus for needed to conduct significance testing on changes in that locus observed between samples."/>
+        <param name="do_normalization" type="select" label="Perform quartile normalization" help="Removes top 25% of genes from FPKM denominator to improve accuracy of differential expression calls for low abundance transcripts.">
+            <option value="No">No</option>
+            <option value="Yes">Yes</option>
+        </param>
        <conditional name="bias_correction">
            <param name="do_bias_correction" type="select" label="Perform Bias Correction" help="Bias detection and correction can significantly improve accuracy of transcript abundance estimates.">
                <option value="Yes">Yes</option>
@@ -103,6 +147,7 @@
                    cuffdiff cuffcompare_out5.gtf cuffdiff_in1.sam cuffdiff_in2.sam 
                -->
                <param name="gtf_input" value="cuffcompare_out5.gtf" ftype="gtf" />
+                <param name="do_groups" value="No" />
                <param name="aligned_reads1" value="cuffdiff_in1.sam" ftype="sam" />
                <param name="aligned_reads2" value="cuffdiff_in2.sam" ftype="sam" />
                <!-- Defaults. -->
@@ -110,6 +155,7 @@
                <param name="min_mapqual" value="0" ftype="sam" />
                <param name="min_alignment_count" value="0" ftype="sam" />
                <param name="do_bias_correction" value="No" />
+                <param name="do_normalization" value="No" />
                <param name="sPaired" value="single" ftype="sam" />
                <!-- 
                    Line diffs are needed because cuffdiff does not produce deterministic output.
@@ -191,6 +237,7 @@ This is a list of implemented Cuffdiff options::
  --FDR FLOAT                    The allowed false discovery rate. The default is 0.05.
  --num-importance-samples INT   Sets the number of importance samples generated for each locus during abundance estimation. Default: 1000
  --max-mle-iterations INT       Sets the number of iterations allowed during maximum likelihood estimation of abundances. Default: 5000
+  -N                             With this option, Cufflinks excludes the contribution of the top 25 percent most highly expressed genes from the number of mapped fragments used in the FPKM denominator. This can improve robustness of differential expression calls for less abundant genes and transcripts.
  
    </help>
 </tool>
@@ -19,7 +19,7 @@ def check_seq_file( dbkey, cached_seqs_pointer_file ):
                seq_path = fields[2].strip()
                break
    return seq_path
-
+	
 def __main__():
    #Parse Command Line
    parser = optparse.OptionParser()
@@ -34,11 +34,14 @@ def __main__():
                                                                                where each end is 50bp, you should set -r to be 200. The default is 45bp.')
    parser.add_option( '-Q', '--min-mapqual', dest='min_mapqual', help='Instructs Cufflinks to ignore alignments with a SAM mapping quality lower than this number. The default is 0.' )
    parser.add_option( '-G', '--GTF', dest='GTF', help='Tells Cufflinks to use the supplied reference annotation to estimate isoform expression. It will not assemble novel transcripts, and the program will ignore alignments not structurally compatible with any reference transcript.' )
+    
+	# Normalization options.
+    parser.add_option( "-N", "--quartile-normalization", dest="do_normalization", action="store_true" )

    # Advanced Options:	
    parser.add_option( '--num-importance-samples', dest='num_importance_samples', help='Sets the number of importance samples generated for each locus during abundance estimation. Default: 1000' )
    parser.add_option( '--max-mle-iterations', dest='max_mle_iterations', help='Sets the number of iterations allowed during maximum likelihood estimation of abundances. Default: 5000' )
-    
+
    # Bias correction options.
    parser.add_option( '-r', dest='do_bias_correction', action="store_true", help='Providing Cufflinks with a multifasta file via this option instructs it to run our new bias detection and correction algorithm which can significantly improve accuracy of transcript abundance estimates.')
    parser.add_option( '', '--dbkey', dest='dbkey', help='The build of the reference dataset' )
@@ -105,6 +108,8 @@ def __main__():
        cmd += ( " --num-importance-samples %i" % int ( options.num_importance_samples ) )
    if options.max_mle_iterations:
        cmd += ( " --max-mle-iterations %i" % int ( options.max_mle_iterations ) )
+    if options.do_normalization:
+        cmd += ( " -N" )
    if options.do_bias_correction:
        cmd += ( " -r %s" % seq_path )
        
@@ -21,6 +21,11 @@
            #if $singlePaired.sPaired == "Yes":
                -m $singlePaired.mean_inner_distance
                -s $singlePaired.inner_distance_std_dev
+                #end if
+
+            ## Normalization?
+            #if str($do_normalization) == "Yes":
+            -N
            #end if
            
            ## Bias correction?
@@ -41,6 +46,10 @@
        <param name="min_isoform_fraction" type="float" value="0.05" min="0" max="1" label="Min Isoform Fraction" help=""/>
        <param name="pre_mrna_fraction" type="float" value="0.05" min="0" max="1" label="Pre MRNA Fraction" help=""/>
        <param name="min_map_quality" type="integer" value="0" min="0" max="255" label="Min SAM Map Quality" help=""/>
+        <param name="do_normalization" type="select" label="Perform quartile normalization" help="Removes top 25% of genes from FPKM denominator to improve accuracy of differential expression calls for low abundance transcripts.">
+            <option value="No">No</option>
+            <option value="Yes">Yes</option>
+        </param>
        <conditional name="reference_annotation">
            <param name="use_ref" type="select" label="Use Reference Annotation">
                <option value="No">No</option>
@@ -102,6 +111,7 @@
            <param name="pre_mrna_fraction" value="0.05"/>
            <param name="min_map_quality" value="0"/>
            <param name="use_ref" value="No"/>
+            <param name="do_normalization" value="No" />
            <param name="do_bias_correction" value="No"/>
            <output name="genes_expression" format="tabular" file="cufflinks_out3.expr"/>
            <output name="transcripts_expression" format="tabular" file="cufflinks_out2.expr"/>
@@ -202,5 +212,6 @@ This is a list of implemented Cufflinks options::
  -j        Some RNA-Seq protocols produce a significant amount of reads that originate from incompletely spliced transcripts, and these reads can confound the assembly of fully spliced mRNAs. Cufflinks uses this parameter to filter out alignments that lie within the intronic intervals implied by the spliced alignments. The minimum depth of coverage in the intronic region covered by the alignment is divided by the number of spliced reads, and if the result is lower than this parameter value, the intronic alignments are ignored. The default is 5%.
  -Q	    Instructs Cufflinks to ignore alignments with a SAM mapping quality lower than this number. The default is 0.
  -G	    Tells Cufflinks to use the supplied reference annotation to estimate isoform expression. It will not assemble novel transcripts, and the program will ignore alignments not structurally compatible with any reference transcript.  
+  -N        With this option, Cufflinks excludes the contribution of the top 25 percent most highly expressed genes from the number of mapped fragments used in the FPKM denominator. This can improve robustness of differential expression calls for less abundant genes and transcripts.
    </help>
 </tool>