| 1 |
#!/usr/bin/perl |
| 2 |
use strict; |
| 3 |
use Getopt::Std; |
| 4 |
use FindBin;use lib $FindBin::Bin; |
| 5 |
|
| 6 |
my $usage = q/Usage: |
| 7 |
rnaseq_cufflinks.pl <libprefix> |
| 8 |
Prepares and runs tophat+cufflinks on a pair of fastq files |
| 9 |
e.g. |
| 10 |
rnaseq_cufflinks.pl RIVA |
| 11 |
|
| 12 |
/; |
| 13 |
|
| 14 |
my %libs=('BJAB'=>['GCB',250], |
| 15 |
'Dogum'=>['GCB', 250], |
| 16 |
'DOHH2'=>['GCB', 250], |
| 17 |
'Gumbus'=>['GCB',250], |
| 18 |
'HT'=>['GCB',250], |
| 19 |
'Karpas422'=>['GCB',250], |
| 20 |
'Ly1'=>['GCB',150], |
| 21 |
'Ly18'=>['GCB',250], |
| 22 |
'Ly19'=>['GCB',250], |
| 23 |
'Ly2'=>['GCB',250], |
| 24 |
'Ly4'=>['GCB',250], |
| 25 |
'Ly7'=>['GCB',250], |
| 26 |
'Ly8'=>['GCB',150], |
| 27 |
'Mieu'=>['GCB',250], |
| 28 |
'Pfieffer'=>['GCB',250], |
| 29 |
'ULA'=>['GCB',250], |
| 30 |
# ABC: |
| 31 |
'HBL1'=>['ABC',150], |
| 32 |
'TMD8'=>['ABC',150], |
| 33 |
'Ly10'=>['ABC',150], |
| 34 |
'Ly3'=>['ABC',250], |
| 35 |
'RIVA'=>['ABC',250], |
| 36 |
'SUDHL2'=>['ABC',150], |
| 37 |
'U2932'=>['ABC',150] |
| 38 |
); |
| 39 |
|
| 40 |
my $tcmd='tophat --solexa1.3-quals --coverage-search --microexon-search --closure-search -p 8 '; |
| 41 |
#my $tophat1=$tcmd.' -r 36 -o RIVA_tophat_untrimmed hg19 RIVA_1.fq RIVA_2.fq > & log_tophat_untrimmed.log &'; |
| 42 |
#my $tophat2=$tcmd.' -r 80 -o RIVA_cut84_tophat hg19 RIVA_1.unmapped_cut84.fq RIVA_2.unmapped_cut84.fq > & log_tophat_cut84.log &'; |
| 43 |
my $gzdir='/fs/szattic-asmg6/NCI_rnaseq/Cell_line_seq'; |
| 44 |
my $workdir='/fs/szannotation/human_rnaseq/NCI'; |
| 45 |
umask 0002; |
| 46 |
chdir($workdir) || die("Error switching to working directory: $workdir\n"); |
| 47 |
|
| 48 |
getopts('o:') || die($usage."\n"); |
| 49 |
my $outfile=$Getopt::Std::opt_o; |
| 50 |
if ($outfile) { |
| 51 |
open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n"); |
| 52 |
select(OUTF); |
| 53 |
} |
| 54 |
# -- |
| 55 |
|
| 56 |
my $fbase=$ARGV[0]; |
| 57 |
my $libd=$libs{$fbase}; |
| 58 |
die("$usage Error: unrecognized sample ($fbase)!\n") unless $libd; |
| 59 |
my ($gz1, $gz2)=($fbase.'_1_sequence.txt.gz',$fbase.'_2_sequence.txt.gz'); |
| 60 |
foreach my $gz ($gz1, $gz2) { |
| 61 |
die("Error: file $gzdir/$gz not found!\n") unless -f "$gzdir/$gz"; |
| 62 |
} |
| 63 |
my ($fq1, $fq2)=($fbase.'_1.fq', $fbase.'_2.fq'); |
| 64 |
my ($fqu1, $fqu2)=($fbase.'_unmapped_1.fq', $fbase.'_unmapped_2.fq'); |
| 65 |
my ($fqustats1, $fqustats2)=($fbase.'_unmapped_1.qstats', $fbase.'_unmapped_2.qstats'); |
| 66 |
my ($fqcut1, $fqcut2)=($fbase.'_trimmed_1.fq', $fbase.'_trimmed_2.fq'); |
| 67 |
if (-f $fq1) { |
| 68 |
mlog("$fq1 already exists, skipping decompression."); |
| 69 |
} |
| 70 |
else { |
| 71 |
mlog("decompressing $gz1 .."); |
| 72 |
system("gzip -cd $gzdir/$gz1 > $fq1"); |
| 73 |
} |
| 74 |
unless (-f "$fq1.cidx") { |
| 75 |
mlog("indexing $fq1.."); |
| 76 |
system("cdbfasta -Q $fq1"); |
| 77 |
system("cdbyank -l $fq1.cidx | sort -u > $fq1.lst"); |
| 78 |
} |
| 79 |
# fq2 |
| 80 |
if (-f $fq2) { |
| 81 |
mlog("$fq2 already exists, skipping decompression."); |
| 82 |
} |
| 83 |
else { |
| 84 |
mlog("decompressing $gz2 .."); |
| 85 |
system("gzip -cd $gzdir/$gz2 > $fq2"); |
| 86 |
} |
| 87 |
unless (-f "$fq2.cidx") { |
| 88 |
mlog("indexing $fq2.."); |
| 89 |
system("cdbfasta -Q $fq2"); |
| 90 |
system("cdbyank -l $fq2.cidx | sort -u > $fq2.lst"); |
| 91 |
} |
| 92 |
unless (-s "$fbase.lst" ) { |
| 93 |
system("sort -o $fbase.lst -m $fq1.lst $fq2.lst"); |
| 94 |
} |
| 95 |
my $rlen=int(`head -2 $fq1 | tail -1 | wc -L`); |
| 96 |
my $frgsize=$libd->[1]; # |
| 97 |
my $innerdist=$frgsize-$rlen*2; |
| 98 |
mlog("read length: $rlen, fragment len: $frgsize (inner dist=$innerdist)"); |
| 99 |
#die("Just for now, exit\n"); |
| 100 |
my $tophat_dir=$fbase.'_tophat_untrimmed'; |
| 101 |
if (-d $tophat_dir) { |
| 102 |
mlog("dir $tophat_dir exists, skipping first tophat run."); |
| 103 |
} |
| 104 |
else { |
| 105 |
mkdir($tophat_dir) || die("Error creating directory $tophat_dir ($!)\n"); |
| 106 |
my $tophatcmd="$tcmd -r $innerdist -o $tophat_dir hg19 $fq1 $fq2 >& $tophat_dir/tophat.log"; |
| 107 |
mlog('Starting tophat #1 run on the untrimmed data:', ' '.$tophatcmd); |
| 108 |
system($tophatcmd); |
| 109 |
mlog('tophat run #1 done.'); |
| 110 |
} |
| 111 |
unless (-s "$tophat_dir/mapped_reads.lst") { |
| 112 |
mlog("generating the list of mapped reads ($tophat_dir/mapped_reads.lst"); |
| 113 |
sysrun("samtools view $tophat_dir/accepted_hits.bam | sam2reads.pl | sort -u > $tophat_dir/mapped_reads.lst"); |
| 114 |
} |
| 115 |
# - generate the fq file with the unmapped reads |
| 116 |
unless (-f $fqu1) { |
| 117 |
mlog("extracting unmapped reads from $fq1"); |
| 118 |
sysrun("comm -13 $tophat_dir/mapped_reads.lst $fq1.lst | cdbyank $fq1.cidx > $fqu1"); |
| 119 |
} |
| 120 |
unless (-f $fqu2) { |
| 121 |
mlog("extracting unmapped reads from $fq2"); |
| 122 |
sysrun("comm -13 $tophat_dir/mapped_reads.lst $fq2.lst | cdbyank $fq2.cidx > $fqu2"); |
| 123 |
} |
| 124 |
#analyse the qv distribution for $fqu1 and $fqu2, find the cutoff point |
| 125 |
unless (-f $fqustats1) { |
| 126 |
mlog("getting quality stats for $fqu1"); |
| 127 |
sysrun("fastx_quality_stats -i $fqu1 -o $fqustats1"); |
| 128 |
} |
| 129 |
|
| 130 |
unless (-f $fqustats2) { |
| 131 |
mlog("getting quality stats for $fqu2"); |
| 132 |
sysrun("fastx_quality_stats -i $fqu2 -o $fqustats2"); |
| 133 |
} |
| 134 |
|
| 135 |
my ($q1cut, $q2cut); |
| 136 |
open(FQ, $fqustats1); |
| 137 |
while (<FQ>) { |
| 138 |
my @t=split(); |
| 139 |
next unless $t[0]>5; |
| 140 |
$q1cut=$t[0]; |
| 141 |
last if $t[5]<10 |
| 142 |
} |
| 143 |
close(FQ); |
| 144 |
|
| 145 |
open(FQ, $fqustats2); |
| 146 |
while (<FQ>) { |
| 147 |
my @t=split(); |
| 148 |
next unless $t[0]>5; |
| 149 |
$q2cut=$t[0]; |
| 150 |
last if $t[5]<10 |
| 151 |
} |
| 152 |
close(FQ); |
| 153 |
|
| 154 |
my $cutlen= ($q1cut<$q2cut ? $q1cut : $q2cut); |
| 155 |
die("Error getting cut length ($q1cut, $q2cut)!\n") if $cutlen<20; |
| 156 |
unless ( -f $fqcut1 ) { |
| 157 |
mlog("trimming $fqu1 to length $cutlen (into $fqcut1)"); |
| 158 |
sysrun("fastx_trimmer -f 1 -l $cutlen -i $fqu1 -o $fqcut1"); |
| 159 |
} |
| 160 |
unless ( -f $fqcut2 ) { |
| 161 |
mlog("trimming $fqu2 to length $cutlen (into $fqcut2)"); |
| 162 |
sysrun("fastx_trimmer -f 1 -l $cutlen -i $fqu2 -o $fqcut2"); |
| 163 |
} |
| 164 |
|
| 165 |
my $tophat_dir2=$fbase.'_tophat_cut'; |
| 166 |
my $cutinnerdist=$frgsize-$cutlen*2; |
| 167 |
mlog("cut read length: $cutlen, fragment len: $frgsize (inner dist=$cutinnerdist)"); |
| 168 |
|
| 169 |
if (-d $tophat_dir2) { |
| 170 |
mlog("dir $tophat_dir2 exists, skipping 2nd tophat run."); |
| 171 |
} |
| 172 |
else { |
| 173 |
mkdir($tophat_dir2) || die("Error creating directory $tophat_dir ($!)\n"); |
| 174 |
my $tophatcmd="$tcmd -r $cutinnerdist -o $tophat_dir2 hg19 $fqcut1 $fqcut2 >& $tophat_dir2/tophat.log"; |
| 175 |
mlog('Starting tophat #2 run on the trimmed data:', ' '.$tophatcmd); |
| 176 |
system($tophatcmd); |
| 177 |
mlog("tophat run #2 done (exit code $?)"); |
| 178 |
} |
| 179 |
#join resulting bam files for cufflinks |
| 180 |
chdir($tophat_dir) || die("Error at chdir into $tophat_dir!\n"); |
| 181 |
unless (-f 'merged_hits.bam') { |
| 182 |
mlog("merging mappings into $tophat_dir/merged_hits.bam"); |
| 183 |
sysrun("samtools merge merged_hits.bam accepted_hits.bam ". |
| 184 |
"../$tophat_dir2/accepted_hits.bam"); |
| 185 |
} |
| 186 |
unless (-d 'cufflinks_out') { |
| 187 |
my $idist=int(($innerdist+$cutinnerdist)/2); |
| 188 |
#my $cufflinks_cmd="cufflinks.new -m $idist -s 34 -p 8 -o cufflinks_out ". |
| 189 |
# ' merged_hits.bam >& log_cufflinks.log'; |
| 190 |
my $cufflinks_cmd="cufflinks -m $frgsize -s 60 -p 8 -o cufflinks_out ". |
| 191 |
' merged_hits.bam >& log_cufflinks.log'; |
| 192 |
mlog("running cufflinks command:\n$cufflinks_cmd"); |
| 193 |
system($cufflinks_cmd); |
| 194 |
mlog("cufflinks done (exit code $?)"); |
| 195 |
} |
| 196 |
# -- |
| 197 |
if ($outfile) { |
| 198 |
select(STDOUT); |
| 199 |
close(OUTF); |
| 200 |
} |
| 201 |
print STDERR "done.\n"; |
| 202 |
#************ Subroutines ************** |
| 203 |
|
| 204 |
sub mlog { |
| 205 |
foreach my $m (@_) { |
| 206 |
print STDERR $m."\n"; |
| 207 |
} |
| 208 |
} |
| 209 |
|
| 210 |
sub sysrun { |
| 211 |
my ($cmd)=@_; |
| 212 |
my $errmsg = `($cmd) 2>&1`; |
| 213 |
my $exitcode=$?; |
| 214 |
if ($exitcode || ($errmsg=~/Error|Segmentation|Fail|Invalid|Cannot|Unable/si)) { |
| 215 |
print STDERR "!Error at:\n$cmd\n"; |
| 216 |
print STDERR " Exit code: $exitcode, message:\n$errmsg\n"; |
| 217 |
#unlink(@todel); |
| 218 |
exit(1); |
| 219 |
} |
| 220 |
} |
