annotate pe_sync_2_files.pl @ 1:b0ab279b5add default tip

Add test cases
author Jim Johnson <jj@umn.edu>
date Wed, 20 Mar 2013 08:26:48 -0500
parents 751f4938cf0d
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
1 #!/usr/bin/perl -w
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
2
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
3 ###############################################
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
4 # pe-sync-2-files.pl # this is a stripped-down version of pe-sync-directory.pl
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
5 # John Garbe
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
6 # February 2012
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
7 #
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
8 # Paired-end read file synchronization script: checks if left and right
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
9 # reads are in the same order.
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
10 # This script checks the two files specificd on the command line.
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
11 # The "quick" option will not report the percent of out-of-sync reads
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
12 # for many failed files, but will run much faster.
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
13 #
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
14 # This script should run fine on linux systems or solaris
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
15 # (which is what the project space fileserver runs)
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
16 #
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
17 # This script can handle Casava 1.8.0 style read IDs, and pre 1.8.0 style ids.
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
18 # Other types of read ID formats will cause a terminal error.
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
19 #
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
20 ###############################################
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
21
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
22 die "Usage: pe-sync-directory.pl [quick] left_fastqfile right_fastqfile\n" unless ($#ARGV > 0);
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
23
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
24 # determine if we're running on solaris in order to change system calls
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
25 $output = `perl -v`;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
26 $solaris = 0;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
27 $solaris = 1 if ($output =~ /solaris/);
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
28 if ($solaris) {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
29 # print STDERR "Running on Solaris\n";
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
30 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
31
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
32 # handle "quick" option
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
33 $quick = 0; # 1 for quickcheck, 0 for thorough check
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
34 if (defined($ARGV[0])) {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
35 if ($ARGV[0] eq "quick") {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
36 $quick = 1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
37 shift @ARGV;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
38 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
39 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
40 print STDERR "QUICK CHECK enabled\n" if ($quick);
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
41
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
42 $filename_r1 = $ARGV[0];
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
43 $filename_r2 = $ARGV[1];
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
44
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
45 # identify read id format type
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
46 if ($solaris) {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
47 $readid = `head -1 $filename_r1`;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
48 } else {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
49 $readid = `head -n 1 $filename_r1`;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
50 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
51 chomp $readid;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
52 if ($readid =~ /^@\S+\/[12]$/) { # @ at the start of the line followed by non-whitespace, a /, a 1 or 2, the end of the line
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
53 $id_type = 1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
54 print STDERR "Casava 1.7 read id style\n"; # TESTING
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
55 } elsif ($readid =~ /^@\S+\W[12]\S+$/) { # @ at the start of the line followed by non-whitspace, a space, a 1 or 2, non-whitespace
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
56 $id_type = 2;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
57 print STDERR "Casava 1.8 read id style\n"; # TESTING
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
58 } else {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
59 print STDERR "Cannot determine read id style\n";
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
60 print STDOUT "Unknwon id style: $readid\n";
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
61 exit 1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
62 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
63
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
64 chomp $filename_r1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
65 chomp $filename_r2;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
66
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
67 if ($quick) {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
68 $failure = &quickcheck($filename_r1, $filename_r2); # quickie check
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
69 print STDERR "FAILED quickcheck\n" if ($failure);
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
70 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
71 # If we aren't doing the quick check, or the quick check passed, do the slow check
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
72 if (!($quick) || !($failure)) {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
73 $failure = &percent($filename_r1, $filename_r2); # slow thorough check
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
74 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
75
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
76 exit;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
77
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
78 ######################## subs #########################
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
79 # do a quick check to see if the files are bad (many bad files will pass this quick check)
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
80 sub quickcheck {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
81 my ($file1, $file2) = @_;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
82 if ($solaris) {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
83 $result1 = `tail -10000l $file1 | head -1`; # grab the line 10,000 lines from the EOF
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
84 $result2 = `tail -10000l $file2 | head -1`;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
85 } else {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
86 $result1 = `tail -n 10000 $file1 | head -n 1`;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
87 $result2 = `tail -n 10000 $file2 | head -n 1`;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
88 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
89 # process the ids depending on id style
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
90 if ($id_type == 1) { # 1 or 2 at end of id
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
91 chomp $result1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
92 chomp $result2;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
93 chop $result1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
94 chop $result2;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
95 } else { # first word is a unique read id, second word is left/right-specific
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
96 ($result1) = split ' ', $result1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
97 ($result2) = split ' ', $result2;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
98 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
99 if ($result1 eq $result2) {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
100 # print "pass\n";
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
101 return 0;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
102 } else {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
103 # print "fail\n";
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
104 return 1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
105 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
106 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
107
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
108
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
109 # compare every read id from each file and determine what percent of reads are
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
110 # out of sync
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
111 sub percent {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
112 my($file1, $file2) = @_;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
113
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
114 open R1, '<', $file1 or die "Couldn't open $file1 for reading:$!";
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
115 open R2, '<', $file2 or die "Couldn't open $file2 for reading:$!";
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
116
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
117 my $readCount = 0;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
118 my $failcount = 0;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
119 while (!eof(R1) and !eof(R2)) {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
120 $id1 = <R1>;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
121 $id2 = <R2>;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
122 next unless ($. % 4 == 1); # check every sequence (4th line in the file)
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
123 # process the ids depending on id style
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
124 if ($id_type == 1) { # 1 or 2 at end of id
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
125 chomp $id1; # newline
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
126 chop $id1; # last char of the id (should be the "1" or "2")
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
127 chomp $id2;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
128 chop $id2;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
129 } else {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
130 ($id1) = split ' ', $id1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
131 ($id2) = split ' ', $id2;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
132 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
133 $failcount++ unless ($id1 eq $id2);
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
134 $readCount++;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
135 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
136 my $percent = $failcount / $readCount;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
137 $prettypercent = $percent * 100; # make it percent
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
138 $prettypercent = sprintf("%.2f", $prettypercent);
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
139 if ($failcount > 0) {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
140 print STDERR "FAILED full check\t$prettypercent% of reads are out-of-sync\n";
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
141 $exitStatus=1;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
142 } else {
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
143 print STDERR "PASSED full check\n";
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
144 $exitStatus=0;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
145 }
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
146
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
147 return $exitStatus;
751f4938cf0d Uploaded
jjohnson
parents:
diff changeset
148 }