\n";
-close(HTML);
-
diff -r a1f084714800 -r 4ca4e9200808 edgeR.xml
--- a/edgeR.xml Thu Dec 26 05:36:07 2013 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,214 +0,0 @@
-
- - Estimates differential gene expression for short read sequence count using methods appropriate for count data
-
- edgeR
- limma
-
-
- edgeR.pl -a $analysis_type.analysis -e $html_file.files_path -f $fdr -h $html_file -o $output
- ## Pairwise comparisons
- #if $analysis_type.analysis == "pw":
- -r $analysis_type.rowsumfilter
- #if $analysis_type.tagwise_disp.twd == "TRUE":
- -u $analysis_type.tagwise_disp.twd_trend
- -t
- #end if
- ## GLM
- #else if $analysis_type.analysis == "glm":
- #if $analysis_type.exp.export_norm == "true":
- -n $norm_exp
- #end if
- -d $analysis_type.disp
- $analysis_type.cont_pw
- #for $fct in $analysis_type.factors:
- factor::${$fct.fact_name}::${$fct.fact}
- #end for
- #for $c in $analysis_type.cont_pred:
- cp::${c.cp_name}::${c.cp}
- #end for
- #for $cnt in $analysis_type.contrasts:
- "cnt::${cnt.add_cont}"
- #end for
- ## LIMMA
- #else
- #if $analysis_type.exp.export_norm == "true":
- -n $norm_exp $analysis_type.exp.log
- #end if
- $analysis_type.cont_pw
- #for $fct in $analysis_type.factors:
- factor::${$fct.fact_name}::${$fct.fact}
- #end for
- #for $c in $analysis_type.cont_pred:
- cp::${c.cp_name}::${c.cp}
- #end for
- #for $cnt in $analysis_type.contrasts:
- "cnt::${cnt.add_cont}"
- #end for
- #end if
- $matrix
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- analysis_type[ "analysis" ] != "pw" and analysis_type[ "exp" ][ "export_norm" ] == "true"
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-Estimates differential gene expression for short read sequence count using methods appropriate for count data.
-If you have paired data you may also want to consider Tophat/Cufflinks.
-Input must be raw count data for each sequence arranged in a rectangular matrix as a tabular file.
-Note - no scaling - please make sure you have untransformed raw counts of reads for each sequence.
-
-Performs digital differential gene expression analysis between groups (eg a treatment and control).
-Biological replicates provide information about experimental variability required for reliable inference.
-
-**What it does not do**
-edgeR_ requires biological replicates.
-Without replicates you can't account for known important experimental sources of variability that the approach implemented here requires.
-
-
-**Input**
-A count matrix containing sequence names as rows and sample specific counts of reads from this sequence as columns.
-The matrix must have 2 header rows, the first indicating the group assignment and the second uniquely identifiying the samples. It must also contain a unique set of (eg Feature) names in the first column.
-
-Example::
-
- # G1:Mut G1:Mut G1:Mut G2:WT G2:WT G2:WT
- #Feature Spl1 Spl2 Spl3 Spl4 Spl5 Spl6
- NM_001001130 97 43 61 34 73 26
- NM_001001144 25 8 9 3 5 5
- NM_001001152 72 45 29 20 31 13
- NM_001001160 0 1 1 1 0 0
- NM_001001177 0 1 0 4 3 3
- NM_001001178 0 2 1 0 4 0
- NM_001001179 0 0 0 0 0 2
- NM_001001180 0 0 0 0 0 2
- NM_001001181 415 319 462 185 391 155
- NM_001001182 1293 945 987 297 938 496
- NM_001001183 5 4 11 7 11 2
- NM_001001184 135 198 178 110 205 64
- NM_001001185 186 1 0 1 1 0
- NM_001001186 75 90 91 34 63 54
- NM_001001187 267 236 170 165 202 51
- NM_001001295 5 2 6 1 7 0
- NM_001001309 1 0 0 1 2 1
- ...
-
-
-Please use the "Count reads in features with htseq-count" tool to generate the count matrix.
-
-**Output**
-
-A tabular file containing relative expression levels, statistical estimates of differential expression probability, R scripts, log, and some helpful diagnostic plots.
-
-.. class:: infomark
-
-**Attribution**
-This tool wraps the edgeR_ Bioconductor package so all calculations and plots are controlled by that code. See edgeR_ for all documentation and appropriate attribution.
-Recommended reference is Mark D. Robinson, Davis J. McCarthy, Gordon K. Smyth, PMCID: PMC2796818
-
-.. class:: infomark
-
-**Attribution**
-When applying the LIMMA (Linear models for RNA-Seq) anlysis the tool also makes use of the limma_ Bioconductor package.
-Recommended reference is Smyth, G. K. (2005). Limma: linear models for microarray data. In: 'Bioinformatics and Computational Biology Solutions using R and Bioconductor'. R. Gentleman, V. Carey, S. Dudoit, R. Irizarry, W. Huber (eds), Springer, New York, pages 397--420.
-
- .. _edgeR: http://www.bioconductor.org/packages/release/bioc/html/edgeR.html
- .. _limma: http://www.bioconductor.org/packages/release/bioc/html/limma.html
-
-
-
-
-
diff -r a1f084714800 -r 4ca4e9200808 htseq.pl
--- a/htseq.pl Thu Dec 26 05:36:07 2013 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,107 +0,0 @@
-#!/usr/bin/perl
-
-use strict;
-use warnings;
-use Getopt::Std;
-use File::Basename;
-$| = 1;
-
-# Grab and set all options
-my %OPTIONS = (a => 0, i => "gene_id", m => "intersection-nonempty", s => "no", t => "exon");
-getopts('a:cg:i:m:o:r:s:t:', \%OPTIONS);
-
-die qq(
-Usage: HTSeq.pl [OPTIONS] Group1=sample1= [Group1=sample2= ... Group2=sampleN= ...]
-
-OPTIONS: -a STR skip all reads with alignment quality lower than the given minimum value (default: $OPTIONS{a})
- -c reduce the matrix by removing any feature with no counts
- -g STR the features file in the GFF/GTF format
- -i STR GFF attribute to be used as feature ID (default: $OPTIONS{i})
- -m STR mode to handle reads overlapping more than one feature. Possible values for are union, intersection-strict and intersection-nonempty (default: $OPTIONS{m})
- -o STR output file name for expression matrix
- -r STR the name of the output report
- -s STR whether the data is from a strand-specific assay (default: $OPTIONS{s})
- -t STR feature type (3rd column in GFF file) to be used, all features of other type are ignored (default, suitable for RNA-Seq and Ensembl GTF files: $OPTIONS{t})
-
-) if(@ARGV == 0);
-
-my $sam_out;
-my @counts;
-my @features;
-my %report;
-my @samplenames;
-my $current_group;
-my @groups;
-my @files;
-my $groupcount = 0;
-my %grouphash;
-
-foreach my $input (@ARGV) {
- my ($group, $sample, $input) = split "::", $input;
- if(! defined $grouphash{$group}) {
- $groupcount++;
- $grouphash{$group} = "G${groupcount}:$group";
- }
- push @groups, $grouphash{$group};
- push @samplenames, $sample;
- push @files, $input;
-}
-
-for(my $index = 0; $index <= $#files; $index++) {
- my $input_file = $files[$index];
- my $sample = $samplenames[$index];
-
- # run htseq
- my @htseq;
- my $COMM;
- my $file_type = `file $input_file`;
- if(grep /text$/, $file_type ) {
- $COMM = "htseq-count -q -m $OPTIONS{m} -s $OPTIONS{s} -a $OPTIONS{a} -t $OPTIONS{t} -i $OPTIONS{i} $input_file $OPTIONS{g}";
- @htseq = `$COMM`;
- } else {
- $COMM = "samtools view $input_file | htseq-count -q -m $OPTIONS{m} -s $OPTIONS{s} -a $OPTIONS{a} -t $OPTIONS{t} -i $OPTIONS{i} - $OPTIONS{g}";
- @htseq = `$COMM`;
- }
-
- my $row = 0;
- $report{$sample} = "Command Used: $COMM\n";
-
- for(my $row = 0; $row <= $#htseq; $row++) {
- # store the report is an hash
- if(grep /^no_feature|^ambiguous|^too_low_aQual|^not_aligned|^alignment_not_unique/, $htseq[$row]) {
- $report{$sample} .= $htseq[$row];
- } else {
- # store the counts in a matrix
- chomp $htseq[$row];
- my ($feature, $value) = split "\t", $htseq[$row];
- $counts[$row][$index] = $value;
- if($input_file eq $files[0]) {
- push @features, $feature;
- }
- }
- }
-}
-
-# print the matrix
-open(MATRIX, ">$OPTIONS{o}") || die "Could Not Create Output File $OPTIONS{o}!\n";
-print MATRIX "#\t".join("\t", @groups)."\n";
-print MATRIX "#Feature\t".join("\t", @samplenames)."\n";
-for(my $row = 0; $row <= $#features; $row++) {
- if(defined $OPTIONS{c}) {
- my $rowsum = 0;
- $rowsum += $_ foreach @{ $counts[$row] };
- if(!$rowsum) {
- next;
- }
- }
- print MATRIX "$features[$row]\t".join("\t", @{ $counts[$row] })."\n";
-}
-close(MATRIX);
-
-# print the report
-open(REPORT, ">$OPTIONS{r}") || die "Could Not Create Output File $OPTIONS{r}!\n";
-print REPORT "$groups[$_]:$samplenames[$_]\n$report{$samplenames[$_]}\n" foreach (0..$#samplenames);
-close(REPORT);
-
-
-
diff -r a1f084714800 -r 4ca4e9200808 htseq.xml
--- a/htseq.xml Thu Dec 26 05:36:07 2013 -0500
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,105 +0,0 @@
-
-
- samtools
- htseq
-
- - Create a digital expression matrix by counting reads in features with htseq-count
-
- htseq.pl -m $MODE -s $STRANDED -a $MINAQUAL -t $FEATURETYPE -i $IDATTR -g $gff_file -o $output -r $report $REDUCE
- ## Inputs.
- #for $group in $group_analysis
- ${group.group}::${group.sample_init}::${group.file_init}
- #for $input_files in $group.input_files:
- ${group.group}::${input_files.sample}::${input_files.file}
- #end for
- #end for
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-Create a digital expression matrix by counting reads in features with htseq-count
-
-For each given file with aligned sequencing reads this tool counts how many reads map to each feature. It then constructs a matrix where the rows represent the features and the columns represent the files.
-
-A feature is here an interval (i.e., a range of positions) on a chromosome or a union of such intervals.
-
-In the case of RNA-Seq, the features are typically genes, where each gene is considered here as the union of all its exons. One may also consider each exon as a feature, e.g., in order to check for alternative splicing. For comparative ChIP-Seq, the features might be binding region from a pre-determined list.
-
-Special care must be taken to decide how to deal with reads that overlap more than one feature. The htseq-count script allows to choose between three modes.
-
-The three overlap resolution modes of htseq-count work as follows. For each position i in the read, a set S(i) is defined as the set of all features overlapping position i. Then, consider the set S, which is (with i running through all position within the read)
-
- * the union of all the sets S(i) for mode union.
- * the intersection of all the sets S(i) for mode intersection-strict.
- * the intersection of all non-empty sets S(i) for mode intersection-nonempty.
-
-If S contains precisely one feature, the read is counted for this feature. If it contains more than one feature, the read is counted as ambiguous (and not counted for any features), and if S is empty, the read is counted as no_feature.
-
-The following figure illustrates the effect of these three modes:
-
-.. image:: http://www-huber.embl.de/users/anders/HTSeq/doc/_images/count_modes.png
-
-The strandedness of the assay may also be set. For stranded=no, a read is considered overlapping with a feature regardless of whether it is mapped to the same or the opposite strand as the feature. For stranded=yes and single-end reads, the read has to be mapped to the same strand as the feature. For paired-end reads, the first read has to be on the same strand and the second read on the opposite strand. For stranded=reverse, these rules are reversed.
-
-.. class:: warningmark
-
-**Important:** The default for strandedness is no. If yes or reverse is selected and your RNA-Seq data has not been made with a strand-specific protocol, this will cause half of the reads to be lost. Hence, make sure to set the option --stranded=no unless you have strand-specific data!
-
-------
-
-**Output**
-
-The script outputs a digital expression matrix containing the counts for each feature by each input Sam/Bam file. It will also generate a report containing special counters, which count reads that were not counted for any feature for various reasons, namely:
-
- * no_feature: reads which could not be assigned to any feature (set S as described above was empty).
- * ambiguous: reads which could have been assigned to more than one feature and hence were not counted for any of these (set S had more than one element).
- * too_low_aQual: reads which were not counted due to the -a option, see below
- * not_aligned: reads in the Sam/Bam file without alignment
- * alignment_not_unique: reads with more than one reported alignment. These reads are recognized from the NH optional SAM field tag. (If the aligner does not set this field, multiply aligned reads will be counted multiple times.)
-
-**Reference**
-
-http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html
-
-
-
-
-
-