perl/etc/justonce/oyk/FGI/OYK.pm

   1 package main;
   2 use strict;
   3 use warnings;
   4 use Data::Dump qw(ddx);
   5
   6 my $SgeQueue = 'fgi.q';
   7 my $SgeProject = 'fgi';
   8 my $SgeJobPrefix = 'nip'.(time() % 9999);
   9
  10 sub Scutadapt($$$$) {
  11         my ($cwd,$cnt,$flst,$outP) = @_;
  12         my $SHcutadapt = <<"END_SH";
  13 #!/bin/bash
  14 #\$ -S /bin/bash
  15 #\$ -q $SgeQueue -P $SgeProject
  16 #\$ -N p0${SgeJobPrefix}cut
  17 #\$ -l vf=600M,num_proc=2
  18 #\$ -binding linear:3
  19 #\$ -cwd -r y
  20 #\$ -v PERL5LIB,PATH,LD_LIBRARY_PATH
  21 #\$ -e /dev/null -o /dev/null
  22
  23 #\$ -t 1-$cnt
  24
  25 cd $cwd
  26
  27 if [ -t 1 && "x\$SKIP" = 'x1' ]; then
  28         while read -u 9 THELINE; do
  29                 read -ra INDAT <<<"\$THELINE"
  30                 if [ -e "\${INDAT[1]}.fq.gz" ]; then
  31                         ln -si \${INDAT[1]}.fq.gz $outP/\${INDAT[0]}.fq.gz
  32                 elif [ -e "\${INDAT[1]}_1.fq.gz" ] && [ -e "\${INDAT[1]}_2.fq.gz" ]; then
  33                         echo "[!]Cannot use SKIP mode for PE, using Read 1 only [\${INDAT[1]}_1.fq.gz] !"
  34                         ln -si \${INDAT[1]}_1.fq.gz $outP/\${INDAT[0]}.fq.gz
  35                 else
  36                         echo "[x]Cannot find .fq.gz file(s) for [\${INDAT[1]}(_[12])?.fq.gz] !"
  37                 fi
  38         done 9<$flst
  39         exit 0
  40 fi
  41
  42 INFILE=`sed -n "\${SGE_TASK_ID}p" $flst`
  43 read -ra INDAT <<<"\$INFILE"
  44
  45 ADAPTERSTR='-g GAACGACATGGCTACGATCCGACTT -a AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA -A AAGTCGGATCGTAGCCATGTCGTTC -G TTGTCTTCCTAAGACCGCTTGGCCTCCGACTT'
  46
  47 if [ -e "\${INDAT[1]}.fq.gz" ]; then
  48         cutadapt -j 2 -m 1 -O 4 -n 2 \$ADAPTERSTR --interleaved -o $outP/\${INDAT[0]}.fq.gz \${INDAT[1]}.fq.gz >$outP/\${INDAT[0]}.log
  49 elif [ -e "\${INDAT[1]}_1.fq.gz" ] && [ -e "\${INDAT[1]}_2.fq.gz" ]; then
  50         cutadapt -j 2 -m 1 -O 4 -n 2 \$ADAPTERSTR --interleaved -o $outP/\${INDAT[0]}.fq.gz \${INDAT[1]}_1.fq.gz \${INDAT[1]}_2.fq.gz >$outP/\${INDAT[0]}.log
  51 else
  52         echo "[x]Cannot find .fq.gz file(s) for [\${INDAT[1]}(_[12])?.fq.gz] !"
  53         exit 1
  54 fi
  55 END_SH
  56         return $SHcutadapt;
  57 }
  58
  59 sub Sbwamem($$$$$$) {
  60         my ($cwd,$cnt,$flst,$outP,$FQprefix,$pRef) = @_;
  61         my $SHbwa = <<"END_SH";
  62 #!/bin/bash
  63 #\$ -S /bin/bash
  64 #\$ -q $SgeQueue -P $SgeProject
  65 #\$ -N p1${SgeJobPrefix}bwa -hold_jid p0${SgeJobPrefix}cut
  66 #\$ -l vf=6G,num_proc=7
  67 #\$ -binding linear:12
  68 #\$ -cwd -r y
  69 #\$ -v PERL5LIB,PATH,LD_LIBRARY_PATH
  70 #\$ -e /dev/null -o /dev/null
  71
  72 #\$ -t 1-$cnt
  73
  74 cd $cwd
  75
  76 INFILE=`sed -n "\${SGE_TASK_ID}p" $flst`
  77 read -ra INDAT <<<"\$INFILE"
  78
  79 bwa mem -t 12 -Y $pRef -R "\@RG\\tID:\${INDAT[0]}\\tSM:\${INDAT[0]}" -p $FQprefix/\${INDAT[0]}.fq.gz 2>$outP/\${INDAT[0]}.log | samtools view -bS - | samtools sort -m 2097152 -T $outP/\${INDAT[0]}.tmp -o $outP/\${INDAT[0]}.bam
  80
  81 samtools index $outP/\${INDAT[0]}.bam
  82
  83 END_SH
  84         return $SHbwa;
  85 }
  86
  87 1;