annotate htseq.pl @ 8:734516a21b52 draft

Uploaded
author fcaramia
date Thu, 13 Sep 2012 00:51:13 -0400
parents ebd59bc6855c
children 6324eefd9e91
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
4
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
1 #!/usr/bin/perl
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
2
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
3 use strict;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
4 use warnings;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
5 use Getopt::Std;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
6 use File::Basename;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
7 $| = 1;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
8
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
9 # Grab and set all options
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
10 my %OPTIONS = (a => 0, i => "gene_id", m => "intersection-nonempty", s => "no", t => "exon");
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
11 getopts('a:cg:i:m:o:r:s:t:', \%OPTIONS);
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
12
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
13 die qq(
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
14 Usage: HTSeq.pl [OPTIONS] Group1=sample1=<SAM/BAM file> [Group1=sample2=<SAM/BAM file> ... Group2=sampleN=<SAM/BAM file> ...]
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
15
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
16 OPTIONS: -a STR skip all reads with alignment quality lower than the given minimum value (default: $OPTIONS{a})
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
17 -c reduce the matrix by removing any feature with no counts
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
18 -g STR the features file in the GFF/GTF format
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
19 -i STR GFF attribute to be used as feature ID (default: $OPTIONS{i})
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
20 -m STR mode to handle reads overlapping more than one feature. Possible values for <mode> are union, intersection-strict and intersection-nonempty (default: $OPTIONS{m})
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
21 -o STR output file name for expression matrix
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
22 -r STR the name of the output report
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
23 -s STR whether the data is from a strand-specific assay (default: $OPTIONS{s})
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
24 -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})
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
25
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
26 ) if(@ARGV == 0);
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
27
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
28 my $sam_out;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
29 my @counts;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
30 my @features;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
31 my %report;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
32 my @samplenames;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
33 my $current_group;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
34 my @groups;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
35 my @files;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
36 my $groupcount = 0;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
37 my %grouphash;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
38
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
39 foreach my $input (@ARGV) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
40 my ($group, $sample, $input) = split "::", $input;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
41 if(! defined $grouphash{$group}) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
42 $groupcount++;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
43 $grouphash{$group} = "G${groupcount}:$group";
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
44 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
45 push @groups, $grouphash{$group};
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
46 push @samplenames, $sample;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
47 push @files, $input;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
48 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
49
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
50 for(my $index = 0; $index <= $#files; $index++) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
51 my $input_file = $files[$index];
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
52 my $sample = $samplenames[$index];
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
53
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
54 # run htseq
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
55 my @htseq;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
56 my $COMM;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
57 my $file_type = `file $input_file`;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
58 if(grep /text$/, $file_type ) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
59 $COMM = "htseq-count -q -m $OPTIONS{m} -s $OPTIONS{s} -a $OPTIONS{a} -t $OPTIONS{t} -i $OPTIONS{i} $input_file $OPTIONS{g}";
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
60 @htseq = `$COMM`;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
61 } else {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
62 $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}";
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
63 @htseq = `$COMM`;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
64 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
65
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
66 my $row = 0;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
67 $report{$sample} = "Command Used: $COMM\n";
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
68
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
69 for(my $row = 0; $row <= $#htseq; $row++) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
70 # store the report is an hash
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
71 if(grep /^no_feature|^ambiguous|^too_low_aQual|^not_aligned|^alignment_not_unique/, $htseq[$row]) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
72 $report{$sample} .= $htseq[$row];
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
73 } else {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
74 # store the counts in a matrix
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
75 chomp $htseq[$row];
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
76 my ($feature, $value) = split "\t", $htseq[$row];
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
77 $counts[$row][$index] = $value;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
78 if($input_file eq $files[0]) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
79 push @features, $feature;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
80 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
81 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
82 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
83 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
84
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
85 # print the matrix
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
86 open(MATRIX, ">$OPTIONS{o}") || die "Could Not Create Output File $OPTIONS{o}!\n";
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
87 print MATRIX "#\t".join("\t", @groups)."\n";
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
88 print MATRIX "#Feature\t".join("\t", @samplenames)."\n";
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
89 for(my $row = 0; $row <= $#features; $row++) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
90 if(defined $OPTIONS{c}) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
91 my $rowsum = 0;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
92 $rowsum += $_ foreach @{ $counts[$row] };
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
93 if(!$rowsum) {
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
94 next;
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
95 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
96 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
97 print MATRIX "$features[$row]\t".join("\t", @{ $counts[$row] })."\n";
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
98 }
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
99 close(MATRIX);
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
100
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
101 # print the report
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
102 open(REPORT, ">$OPTIONS{r}") || die "Could Not Create Output File $OPTIONS{r}!\n";
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
103 print REPORT "$groups[$_]:$samplenames[$_]\n$report{$samplenames[$_]}\n" foreach (0..$#samplenames);
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
104 close(REPORT);
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
105
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
106
ebd59bc6855c Uploaded
fcaramia
parents:
diff changeset
107