Davide Cittaro avatar Davide Cittaro committed 18efeb1 Merge

added a simple vcf annotator

Comments (0)

Files changed (5)

+#!/bin/bash
+# simple script to compute alll bwa alignment in a dir.
+run=false
+test "x$1" = "x--run" && run=true
+
+# from original makefile 
+
+BWA=/usr/local/cluster/bin/bwa
+SAMTOOLS=/usr/local/cluster/bin/samtools
+BWAOPT_ALN=
+BWAOPT_PE=
+BWAOPT_SE=
+PICMERGE=/usr/local/cluster/bin/MergeSamFiles.jar
+PICOPTS=
+VALIDATION_STRINGENCY=SILENT
+CREATE_INDEX=true
+MSD=true
+ASSUME_SORTED=true
+
+#
+# to easily change where to write temporary files scratch: 
+#  default: present directory..
+
+
+# to replace with PBS_WKDIR
+cd $PWD
+
+fcid=`tail -1 SampleSheet.csv | cut -d',' -f1 `
+reference=`tail -1 SampleSheet.csv | cut -d',' -f4`
+name=`tail -1 SampleSheet.csv | cut -d',' -f3`
+rindex=`tail -1 SampleSheet.csv | cut -d',' -f5`
+ref_genome=/lustre1/genomes/$reference/bwa/$reference
+
+echo "$fcid"
+echo "$reference"
+
+JOBLIST=
+BAMFILES=
+
+experiment_name=$name
+
+LOCAL_SCRATCH=/lustre2/scratch/${experiment_name}
+mkdir $LOCAL_SCRATCH
+
+echo $experiment_name
+#------------------------------------------------------------------------#
+
+for file in *$experiment_name*R1*.fastq.gz
+do
+	# foreach combination of file write and submit the job
+
+	# identify rindex lane and index
+
+	echo $file
+	lane=`echo $file | rev | cut -d'_' -f 3 | rev`
+	index=`echo $file  | rev| cut -d'_' -f 1 | rev`
+	indx=` echo $index| cut -d'.' -f 1 `
+#	echo $lane, $index ,$indx
+
+	# this is the name:
+
+	R1=$experiment_name"_"$rindex"_"$lane"_R1_"$indx".fastq.gz"
+	R2=$experiment_name"_"$rindex"_"$lane"_R2_"$indx".fastq.gz"
+	R2_final=$experiment_name"_"$rindex"_"$lane"_R2_"$indx
+	R1_final=$experiment_name"_"$rindex"_"$lane"_R1_"$indx
+	R_final=$experiment_name"_"$rindex"_"$lane"_"$indx
+
+#        job_nameR1="gz_"${R1_final:(-11)} 
+#        job_nameR2="gz_"${R2_final:(-11)} 
+#        job_nameR="sampe_"${R_final:(-8)} 
+        job_nameR1="a_${R1_final}"
+        job_nameR2="a_"${R2_final} 
+        job_nameR="s_"${R_final} 
+	#------------------------------------------------------------------------#
+
+	# write the first job script
+	SCRIPT=job.$R1_final
+	cat <<__EOF__> $SCRIPT
+#PBS -l select=1:ncpus=2:mem=16g
+#PBS -N ${job_nameR1:0:15}
+#PBS -M cittaro.davide@hsr.it
+#PBS -P ${experiment_name}
+#PBS -m a
+
+cd $PWD
+
+/usr/bin/time $BWA aln -t 2 $BWAOPT_ALN $ref_genome $R1 >$LOCAL_SCRATCH/tmp.$R1_final
+__EOF__
+
+	if $run
+	then
+		job_id1="$(qsub job.$R1_final)"
+		test "x$job_id1" = "x" && { echo >&2 "*** error while submitting job $SCRIPT" ; exit 1 ; }
+		qstat  $job_id1 || { echo >&2 "*** couldn't check for job $job_id1 (R1)" ; exit 1 ; }
+	fi
+
+	#------------------------------------------------------------------------#
+
+	# write the second job script
+
+	SCRIPT=job.$R2_final
+	cat <<__EOF__> $SCRIPT
+#PBS -l select=1:ncpus=2:mem=16g
+#PBS -N ${job_nameR2:0:15}
+#PBS -M cittaro.davide@hsr.it
+#PBS -P ${experiment_name}
+#PBS -m a
+
+cd $PWD
+
+/usr/bin/time $BWA aln -t 2 $BWAOPT_ALN $ref_genome $R2 >$LOCAL_SCRATCH/tmp.$R2_final
+__EOF__
+
+	if $run
+	then
+		job_id2="$(qsub job.$R2_final)"
+		test "x$job_id2" = "x" && { echo >&2 "*** error while submitting job $SCRIPT" ; exit 1 ; }
+		qstat   $job_id2 || { echo >&2 "*** couldn't check for job $job_id2 (R2)" ; exit 1 ; }
+	fi
+
+	#------------------------------------------------------------------------#
+
+	# create the third job: runs only if previous two successful.
+
+	SCRIPT=job.final.$R_final
+	cat <<__EOF__> $SCRIPT
+#PBS -l select=1:ncpus=2:mem=24g
+#PBS -W depend=afterok:$job_id1:$job_id2
+#PBS -N ${job_nameR:0:15}
+#PBS -M cittaro.davide@hsr.it
+#PBS -P ${experiment_name}
+#PBS -m a
+
+cd $PWD
+
+/usr/bin/time $BWA sampe $BWAOPT_PE -r "@RG\tID:$experiment_name"_"$lane\tPL:illumina\tPU:$fcid\tLB:$experiment_name\tSM:$experiment_name\tCN:CTGB" $ref_genome $LOCAL_SCRATCH/tmp.$R1_final $LOCAL_SCRATCH/tmp.$R2_final  $R1 $R2 >$LOCAL_SCRATCH/tmp.sampe.$R_final
+$SAMTOOLS view -Su $LOCAL_SCRATCH/tmp.sampe.$R_final | $SAMTOOLS sort - $R_final
+__EOF__
+
+	if $run
+	then
+		id_final="$(qsub $SCRIPT)"
+		test "x$id_final" = "x" && { echo >&2 "*** error while submitting $SCRIPT" ; exit 1 ; }
+		qstat   $id_final || { echo >&2 "*** couldn't check for job $id_final (final)" ; exit 1 ; }
+	fi
+
+	test "x$id_final" != "x" && JOBLIST=$JOBLIST:$id_final
+        BAMFILES=$BAMFILES" "I=$R_final".bam"
+
+done
+
+#------------------------------------------------------------------------#
+
+# 
+
+
+JOBLIST=${JOBLIST#:}
+BAMFILES=${BAMFILES#" "}
+
+job_nameC="combin"$experiment_name 
+SCRIPT=job.combine.$R_final
+
+cat <<__EOF__> $SCRIPT
+#PBS -l select=1:ncpus=8:mem=48g
+#PBS -W depend=afterok:$JOBLIST
+#PBS -N ${job_nameC:0:15}
+#PBS -M cittaro.davide@hsr.it
+#PBS -P ${experiment_name}
+#PBS -m a
+#PBS -m e
+
+cd $PWD
+
+/usr/bin/time  java -jar $PICMERGE $BAMFILES \
+        O=$experiment_name.bam \
+        CREATE_INDEX=true \
+        MSD=true \
+        VALIDATION_STRINGENCY=SILENT \
+        ASSUME_SORTED=true \
+        USE_THREADING=true
+
+/usr/local/cluster/bin/samtools flagstat $experiment_name.bam
+/bin/rm -f $LOCAL_SCRATCH/tmp.$R1_final
+/bin/rm -f $LOCAL_SCRATCH/tmp.$R2_final
+/bin/rm -f $LOCAL_SCRATCH/tmp.sampe.$R_final
+
+/bin/mv $experiment_name.bam $experiment_name.bam.lock
+/bin/mv $experiment_name.bai $experiment_name.bai.lock
+
+/bin/rm -f *.bam *.bai
+
+/bin/mv $experiment_name.bam.lock $experiment_name.bam 
+/bin/mv $experiment_name.bai.lock $experiment_name.bai
+
+rm -fr $LOCAL_SCRATCH
+__EOF__
+
+if $run
+then
+	id_combine="$(qsub $SCRIPT)"
+	test "x$id_combine" = "x" && { echo >&2 "*** error while submitting $SCRIPT" ; exit 1 ; }
+	qstat  $id_combine || { echo >&2 "*** couldn't check for job $id_combine (combine)" ; exit 1 ; }
+fi
+
+#------------------------------------------------------------------------#
+
+exit

estimate_distribution.py

   return (r, p)
 
 
-def estimate_distribution(bam_file, bed_file = None, keep_zero = True, upper_limit = 0):
+def estimate_distribution(bam_file, bed_file = None, keep_zero = True, upper_limit = 0, plot_distr = True):
 
   actual_coverage = np.zeros(_max_coverage, dtype=np.double)	# will store the actual coverage
   whole_genome = False											# A flag to differentiate whole genome or targeted regions
   
   # plot actual coverage
   actual_coverage = actual_coverage / np.sum(actual_coverage)
-  plt.plot(actual_coverage, label="Actual Coverage", color='black', linestyle='dashed', marker='o')
+  plt.plot(actual_coverage, label="Actual Coverage", color='black', linestyle=' ', marker='o')
   
   # plot distributions
   stdev = np.sqrt(variance)
   poisson = SS.distributions.poisson(mean)
   nbinom = SS.distributions.nbinom(nb_r, nb_p)
   binom = SS.distributions.binom(b_n, b_p)
-  
+
   yp = poisson.pmf(x)
   yb = binom.pmf(x)
   ynb = nbinom.pmf(x)
   
   ymax = max(yp.max(), yb.max(), ynb.max(), actual_coverage.max())	* 1.1
   
-  plt.plot(x, yp, label="Poisson", color='blue', linestyle='dashed', marker='o')
-  plt.plot(x, yb, label="Binomial", color='green', linestyle='dashed', marker='o')
-  plt.plot(x, ynb, label="Negative Binomial", color='red', linestyle='dashed', marker='o')
+  if plot_distr:
+    plt.plot(x, yp, label="Poisson", color='blue', linestyle=' ', marker='o')
+    plt.plot(x, yb, label="Binomial", color='green', linestyle=' ', marker='o')
+    plt.plot(x, ynb, label="Negative Binomial", color='red', linestyle=' ', marker='o')
   plt.title("Distributions for %s" % (os.path.basename(bam_file)))
   plt.ylim(0, ymax)
   plt.xlim(-1, xmax)
   plt.grid()
   plt.rcParams['legend.fontsize'] = 'small'
   plt.legend(loc=0, ncol=2, scatterpoints=1 , numpoints=1)
-  plt.xlabel('Coverage')
-  plt.ylabel('Probability')
   plt.savefig(os.path.basename(bam_file.replace(".bam", ".png")))
   
   # close fh
   bam_file = None
   keep_zero = True
   upper_limit = 0
-  
+  plot_distr = True
+ 
   # read options
-  options, arguments = getopt.getopt(sys.argv[1:], 'b:vhf:zl:')
+  options, arguments = getopt.getopt(sys.argv[1:], 'b:vhf:zl:n')
   for o, a in options:
     if o == '-b':
       regions_file = os.path.expanduser(a)
       keep_zero = False  
     if o == '-l':
       upper_limit = int(a)  
+    if  o == '-n':
+      plot_distr = False
 
   # some sanity checks:
   if not bam_file or not os.path.exists(bam_file):
     sys.stderr.write("%s\n" % errors[2])
     sys.exit(1)
 
-  estimate_distribution(bam_file, regions_file, keep_zero, upper_limit)    
+  estimate_distribution(bam_file, regions_file, keep_zero, upper_limit, plot_distr)    
 
 
 
 if __name__ == '__main__':
   process_arguments()
 
+
+#!/usr/bin/env python2.7
+
+
+import bx.intervals
+import argparse
+import sys
+import vcf
+
+def terminate(msg=None, exit_status=1):
+  sys.stderr.write("%s\n" % msg)
+  sys.exit(exit_status)
+
+def sniff(s):
+  try:
+    int(s)
+    return 'Integer'
+  except ValueError:
+    try:
+      float(s)
+      return 'Float'
+    except ValueError:
+      return 'String'
+
+def cast(s, t):
+  if t == 'Integer':
+    return int(s)
+  elif t == 'Float':
+    return float(s)
+  elif t == 'Flag':
+    return None  
+  else:
+    return s  
+
+def read_gene_list(fh):
+  
+  # Read gene list, which is a tab separated file with gene in the first column.
+  # If one column is present, returns the genes, if multiple columns are present, returns the gene list 
+  # along with the properties. An header should be present, column names will be the properties
+  
+  mark_names = []
+  gene_list = {}
+  
+  # read first line
+  header = fh.readline()
+  mark_names = header.lstrip('#').strip().split()
+  nf = len(mark_names)
+  if nf > 1 and not header.startswith('#'):
+    terminate("Multiple columns in gene list, but no header provided!")
+  elif nf == 1:
+    gene_list[mark_names[0]] = None
+  
+  for line in fh:
+    fields = line.strip().split()
+    if nf == 1:
+      gene_list[fields[0]] = None
+    else:
+      gene_list[fields[0]] = dict(zip(mark_names[1:], fields[1:]))      
+  return gene_list
+
+  
+def read_refgene(fh, gene_list, prefix="MARK", use_gene_symbol = True):
+  
+  # Read sql table (refGene) from UCSC. This should be quite standard, with 
+  # at least cdsStart, cdsEnd, name2 and chrom in the proper positions...
+   
+  # bin	name	chrom	strand	txStart	txEnd	cdsStart	cdsEnd	exonCount	exonStarts	exonEnds	score	name2	cdsStartStat	cdsEndStat	exonFrames
+
+  tree_dict = {}
+  for line in (fh):
+    if line.startswith('#'): continue
+    fields = line.split()
+    chrom = fields[2]
+    cdsStart = int(fields[6])
+    cdsEnd = int(fields[7])
+    if use_gene_symbol:
+      gene_name = fields[12]
+    else:
+      gene_name = fields[1]
+    if not tree_dict.has_key(chrom):
+      tree_dict[chrom] = bx.intervals.IntervalTree()
+    
+    if gene_list.has_key(gene_name):
+      properties = gene_list[gene_name]
+      if not properties:
+        properties = {}
+      tree_dict['__properties__'] = properties.keys()
+      properties['gene'] = gene_name
+      tree_dict[chrom].add(cdsStart, cdsEnd, properties)
+
+  return tree_dict
+
+def annotate_vcf(fh, gene_props, prefix):
+  
+  vcf_file = vcf.Reader(fh)
+  
+  properties = [prefix + x for x in gene_props['__properties__']]
+  
+  if len(properties) == 0:
+    id = prefix
+    type = 'Flag'
+    num = 1
+    desc = 'Value added by markvcf'
+    vcf_file.infos[id] = vcf.parser._Info(id, num, type, desc)
+  else:
+    for id in properties:
+      type = sniff(id)
+      num = '.'
+      desc = 'Value added by markvcf'
+      vcf_file.infos[id] = vcf.parser._Info(id, num, type, desc)
+  
+  vcf_out = vcf.Writer(sys.stdout, vcf_file, lineterminator='\n')
+  # we now traverse the tree. As soon as we get into an overlapping feature
+  # we update the record entry and the INFO in the header
+  for record in vcf_file:
+    chrom = record.CHROM
+    if not gene_props[chrom]:
+      # no need to process, just write
+      vcf_out.write_record(record)
+    else:
+      # get the intersection, it is an array of dictionaries
+      matched_genes = gene_props[chrom].find(record.POS, record.POS)  
+      if len(matched_genes) == 0:
+        vcf_out.write_record(record)
+      else:
+        # ok, this is tough!
+        # we don't want to add redundancies... I didn't removed from intervaltree
+        # as, well, it is easier to manage here :-)
+        added = set([]) # keeps track of redundacies for this SNP
+        for match in matched_genes:
+          gene = match['gene']
+          if gene in added: continue # skip. In principle annotation file should contain one entry per gene...
+          for (p,v) in match.items():
+            # here p is the property name, v its value
+            id = prefix + p # convert to match above...
+            try:
+              mark_value = cast(v, vcf_file.infos[id].type) # this converts back to the proper data type
+              try:
+                # since variants can overlap more than one gene (yes...) they should be handled as arrays
+                record.INFO[id].append(mark_value)
+              except KeyError:
+                record.INFO[id] = [mark_value]  
+            except:
+              # this happens when we only have one gene, no columns in the annotation file
+              record.INFO[prefix] = True
+            added.add(gene)
+          vcf_out.write_record(record)     
+
+def mark_main():
+
+  # parse options
+  option_parser = argparse.ArgumentParser(
+  description="Mark VCF files for gene properties",
+  prog="markvcf",
+  epilog="For any question, write to cittaro.davide@hsr.it")
+  option_parser.add_argument("--version", action="version", version="%(prog)s 0.1")
+  option_parser.add_argument("-v", "--vcf", help="VCF file to annotate", action='store', type=file, required=True)
+  option_parser.add_argument("-g","--genelist", help="File containing gene list", action='store', type=file, required=True)
+  option_parser.add_argument("-r","--refgene", help="RefGene SQL table from UCSC", action='store', type=file, required=True)
+  option_parser.add_argument("-n","--name", help="Annotation name prefix", action='store', default="MARK")
+  options = option_parser.parse_args()
+  
+  # read the genelist and retur
+  gene_list = read_gene_list(options.genelist)
+  
+  # read sql file and create the interval tree with gene properties
+  gene_properties = read_refgene(options.refgene, gene_list)
+  
+  # read the vcf and, at the same time, output annotated entries
+  # I really would like to use PyVCF interface but at time of writing it had a bug that prevents 
+  # reading of my test vcf files...
+  
+  annotate_vcf(options.vcf, gene_properties, options.name)
+  
+  # that's it!
+
+if __name__ == "__main__": mark_main()
+

soapsplice_submit.sh

+#!/bin/bash
+# simple script to compute alll bwa alignment in a dir.
+run=false
+test "x$1" = "x--run" && run=true
+
+# from original makefile 
+
+SSPLICE=/lustre1/tools/bin/soapsplice
+SAMTOOLS=/usr/local/cluster/bin/samtools
+SSPLICEOPT_ALN="-f 2 -q 1 -j 0"
+PICMERGE=/usr/local/cluster/bin/MergeSamFiles.jar
+PICOPTS=
+VALIDATION_STRINGENCY=SILENT
+CREATE_INDEX=true
+MSD=true
+ASSUME_SORTED=true
+
+#
+# to easily change where to write temporary files scratch: 
+#  default: present directory..
+
+
+# to replace with PBS_WKDIR
+cd $PWD
+
+fcid=`tail -1 SampleSheet.csv | cut -d',' -f1 `
+reference=`tail -1 SampleSheet.csv | cut -d',' -f4`
+name=`tail -1 SampleSheet.csv | cut -d',' -f3`
+rindex=`tail -1 SampleSheet.csv | cut -d',' -f5`
+ref_genome=/lustre1/genomes/$reference/SOAPsplice/${reference}.index
+faidx=/lustre1/genomes/$reference/fa/${reference}.fa.fai
+
+echo "$fcid"
+echo "$reference"
+
+JOBLIST=
+BAMFILES=
+
+experiment_name=$name
+
+LOCAL_SCRATCH=/lustre2/scratch/${experiment_name}
+mkdir $LOCAL_SCRATCH
+
+#set lfs 
+
+lfs setstripe -c -1 -i -1 -s 2M ${LOCAL_SCRATCH}
+
+echo $experiment_name
+#------------------------------------------------------------------------#
+
+for file in *$experiment_name*R1*.fastq.gz
+do
+	# foreach combination of file write and submit the job
+
+	# identify rindex lane and index
+
+	echo $file
+	lane=`echo $file | rev | cut -d'_' -f 3 | rev`
+	index=`echo $file  | rev| cut -d'_' -f 1 | rev`
+	indx=` echo $index| cut -d'.' -f 1 `
+#	echo $lane, $index ,$indx
+
+	# this is the name:
+
+	R1=$experiment_name"_"$rindex"_"$lane"_R1_"$indx".fastq.gz"
+	R2=$experiment_name"_"$rindex"_"$lane"_R2_"$indx".fastq.gz"
+	R2_final=$experiment_name"_"$rindex"_"$lane"_R2_"$indx
+	R1_final=$experiment_name"_"$rindex"_"$lane"_R1_"$indx
+	R_final=$experiment_name"_"$rindex"_"$lane"_"$indx
+
+#        job_nameR1="gz_"${R1_final:(-11)} 
+#        job_nameR2="gz_"${R2_final:(-11)} 
+#        job_nameR="sampe_"${R_final:(-8)} 
+        job_nameR="a_"${R_final} 
+	#------------------------------------------------------------------------#
+
+	# write the first job script
+	SCRIPT=job.align.$R_final
+	cat <<__EOF__> $SCRIPT
+#PBS -l select=1:ncpus=4:mem=32g
+#PBS -N a_${R_final:0:13}
+#PBS -M cittaro.davide@hsr.it
+#PBS -P ${experiment_name}
+#PBS -m a
+
+cd $PWD
+
+#create header
+awk '{OFS="\t"; print "@SQ","SN:"\$1,"LN:"\$2}' $faidx > ${LOCAL_SCRATCH}/header.$R_final
+echo -e "@RG\tID:$experiment_name"_"$lane\tPL:illumina\tPU:$fcid\tLB:$experiment_name\tSM:$experiment_name\tCN:CTGB" >> ${LOCAL_SCRATCH}/header.$R_final
+sversion=( \`$SSPLICE | head -n1\` )
+echo -e "@PG\tID:soapsplice\tPN:soapsplice\tVN:\${sversion[2]}" >> ${LOCAL_SCRATCH}/header.$R_final
+
+$SSPLICE -d $ref_genome -1 $R1 -2 $R2 -o ${LOCAL_SCRATCH}/${R_final} -p 4 ${SSPLICEOPT_ALN}
+__EOF__
+
+	if $run
+	then
+		job_id1="$(qsub ${SCRIPT})"
+		test "x$job_id1" = "x" && { echo >&2 "*** error while submitting job $SCRIPT" ; exit 1 ; }
+		qstat  $job_id1 || { echo >&2 "*** couldn't check for job $job_id1 (R1)" ; exit 1 ; }
+	fi
+
+	# create the third job: runs only if previous two successful.
+
+	SCRIPT=job.final.$R_final
+	cat <<__EOF__> $SCRIPT
+#PBS -l select=1:ncpus=2:mem=24g
+#PBS -W depend=afterok:$job_id1
+#PBS -N s_${R_final:0:13}
+#PBS -M cittaro.davide@hsr.it
+#PBS -P ${experiment_name}
+#PBS -m a
+
+cd $PWD
+
+cat ${LOCAL_SCRATCH}/header.$R_final ${LOCAL_SCRATCH}/${R_final}.sam | $SAMTOOLS view -Su - | $SAMTOOLS sort - $R_final
+__EOF__
+
+	if $run
+	then
+		id_final="$(qsub $SCRIPT)"
+		test "x$id_final" = "x" && { echo >&2 "*** error while submitting $SCRIPT" ; exit 1 ; }
+		qstat   $id_final || { echo >&2 "*** couldn't check for job $id_final (final)" ; exit 1 ; }
+	fi
+
+	test "x$id_final" != "x" && JOBLIST=$JOBLIST:$id_final
+        BAMFILES=$BAMFILES" "I=$R_final".bam"
+
+done
+
+#------------------------------------------------------------------------#
+
+# 
+
+
+JOBLIST=${JOBLIST#:}
+BAMFILES=${BAMFILES#" "}
+
+job_nameC="combin"$experiment_name 
+SCRIPT=job.combine.$R_final
+
+cat <<__EOF__> $SCRIPT
+#PBS -l select=1:ncpus=8:mem=48g
+#PBS -W depend=afterok:$JOBLIST
+#PBS -N ${job_nameC:0:15}
+#PBS -M cittaro.davide@hsr.it
+#PBS -P ${experiment_name}
+#PBS -m a
+#PBS -m e
+
+cd $PWD
+
+/usr/bin/time  java -jar $PICMERGE $BAMFILES \
+        O=$experiment_name.bam \
+        CREATE_INDEX=true \
+        MSD=true \
+        VALIDATION_STRINGENCY=SILENT \
+        ASSUME_SORTED=true \
+        USE_THREADING=true
+
+
+# combine junction usage
+cat ${LOCAL_SCRATCH}/*.junc > ${experiment_name}.junc
+
+/usr/local/cluster/bin/samtools flagstat $experiment_name.bam
+/bin/rm -f ${LOCAL_SCRATCH}/${R_final}.sam
+/bin/rm -f ${LOCAL_SCRATCH}/header.$R_final 
+/bin/rm -f ${LOCAL_SCRATCH}/*.junc
+
+/bin/mv $experiment_name.bam $experiment_name.bam.lock
+/bin/mv $experiment_name.bai $experiment_name.bai.lock
+
+/bin/rm -f *.bam *.bai
+
+/bin/mv $experiment_name.bam.lock $experiment_name.bam 
+/bin/mv $experiment_name.bai.lock $experiment_name.bai
+
+rm -fr $LOCAL_SCRATCH
+__EOF__
+
+if $run
+then
+	id_combine="$(qsub $SCRIPT)"
+	test "x$id_combine" = "x" && { echo >&2 "*** error while submitting $SCRIPT" ; exit 1 ; }
+	qstat  $id_combine || { echo >&2 "*** couldn't check for job $id_combine (combine)" ; exit 1 ; }
+fi
+
+#------------------------------------------------------------------------#
+
+exit
+import sys
+import os
+
+fh = open(sys.argv[1])
+
+content = fh.readlines()
+
+sep = '   +------------------------------------+------------+------------+'
+
+print ".. table:: **Stats for %s**" % os.path.basename(sys.argv[1])
+print
+print sep
+print '   |                                    |   number   |  %         |'
+print sep
+for line in content[5:17]:
+  f = line.split(':')
+  v = f[1].split()
+  if len(v) == 1:
+    v.append('')
+  print '   |%-36s|%-12s|%-12s|' % (f[0], v[0], v[1])
+  print sep
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.