Merge pull request #33 from mzhibo/master

Very new to coding and this is my first "big" project in modifying a script. So welcome to make it better!

This PBS version is modified mainly based on the LSF version, but also took other versions as reference. Main changed to create this version is build the job dependencies to guarantee the sequential steps of each job in the original juicer.sh script. #PBS -W depend=afterok:jobID headerline is used.

The launch stats steps and post processing steps are moved into two separate scripts to avoid the variable value problems (expansion) in multiple level of heredocs. These two scripts will be called from the main juicer.sh at the appropriate step.

As a summary, the follwing scripts have been modified:

juicer.sh

split_rmdups.awk

mega.sh

The following two scripts has been added (extracted from the juicer.sh)

launch_stats.sh

postprocessing.sh

All steps has been tested till the ARROWHEAD step.

*HiCCUP step has not been tested in this version because of some un-resolved error in loading Jcuda native libraries.

One important notice for successful run of this PBS version:

***The $groupname variable has to be maximally 7 characters long*****

This length may change according to your cluster setup. See below for the reason.

The jobs' dependency of this PBS version juicer is based on jobID of each step. And the jobID of each job is obtained through the qstat output based on the specific job Name, jobID_stepx= $(qstat |grep <specific job name string> |cut -c 1-7).  Unfortunately, PBS can not use job Name in building job dependency. The job Name column in qstat output has a maximum length of 16 characters by default. When the job Name exceeds this max length, the beginning part of the job Name string will be replaced by “…” and then followed by the last 13 characters (see examples below). If the job name is not fully displayed, the jobID value may be null due to the failure at grep step.  Each job Name contains the $groupname, which gives each run specific timestamp and avoid's disruption from other runs of juicer from the same user. A potential cause of this problem is the ${groupname} value being too long. I have added comments in the juicer.sh script for setting this variable. 

Below is a example of job names being too long:

$qstat

Job ID                    Name             User            Time Use S Queue

------------------------- ---------------- --------------- -------- - -----

2294690.pbs                STDIN            mzhibo          00:00:15 R batch          

2294778.pbs                make4Cbedgraph   mzhibo                 0 R batch          

2294784.pbs                ...dgraph1111111 mzhibo                 0 Q batch          

$

Zhibo

/*

 * The MIT License (MIT)

 *

 * Copyright (c) 2011-2015 Broad Institute, Aiden Lab

 *

 * Permission is hereby granted, free of charge, to any person obtaining a copy

 * of this software and associated documentation files (the "Software"), to deal

 * in the Software without restriction, including without limitation the rights

 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

 * copies of the Software, and to permit persons to whom the Software is

 * furnished to do so, subject to the following conditions:

 *

 * The above copyright notice and this permission notice shall be included in

 * all copies or substantial portions of the Software.

 *

 *  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

 *  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

 *  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

 *  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

 *  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,

 *  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN

 *  THE SOFTWARE.

 *

 * Based on PicardTools 

 * Juicer version 1.5

 */

import java.text.NumberFormat;

import java.text.ParseException;

import java.util.Locale;

import java.io.FileReader;

import java.io.BufferedReader;

import java.io.IOException;

import java.io.File;

class LibraryComplexity {

  public static void main(String[] args) {

    if (args.length != 2 && args.length != 3 && args.length != 1) {

      System.out.println("Usage: java LibraryComplexity <directory> <output file>");

      System.out.println("     : java LibraryComplexity <unique> <pcr> <opt>");

      System.exit(0);

    }

    NumberFormat nf = NumberFormat.getInstance(Locale.US);

    BufferedReader reader = null;

    long readPairs = 0;

    long uniqueReadPairs = 0;

    long opticalDups = 0;

    long totalReadPairs = 0;

    if (args.length == 2 || args.length==1) {

      try {

        File f = new File(args[0] + "/opt_dups.txt");

        if (f.exists()) {

          reader = new BufferedReader(new FileReader(f));

          while (reader.readLine() != null) opticalDups++;

          reader.close();

        }

        f = new File(args[0] + "/merged_nodups.txt");

        if (f.exists()) {

          reader = new BufferedReader(new FileReader(f));

          while (reader.readLine() != null) uniqueReadPairs++;

          reader.close();				

        }

        f = new File(args[0] + "/dups.txt");

        if (f.exists()) {

          reader = new BufferedReader(new FileReader(args[0] + "/dups.txt"));

          while (reader.readLine() != null) readPairs++;

          reader.close();				

          readPairs += uniqueReadPairs; 

        }

        String fname = "inter.txt";

        if (args.length == 2) fname = args[1];

        f = new File(args[0] + "/" + fname);

        if (f.exists()) {

          reader = new BufferedReader(new FileReader(args[0] + "/" + fname));

          String line = reader.readLine();

          boolean done = false;

          while (line != null && !done) {

            if (line.contains("Sequenced Read")) {

              String[] parts = line.split(":");

              try {

                totalReadPairs = nf.parse(parts[1].trim()).longValue();

              }

              catch (ParseException e) {

                totalReadPairs = 0;

              }

              done = true;

            }

            line = reader.readLine(); 

          }

          reader.close();

        }

      } catch (IOException error) {

        System.err.println("Problem counting lines in merged_nodups and dups");

        System.exit(1);

      }

    }

    else {

      try {

        uniqueReadPairs = Integer.valueOf(args[0]);

        readPairs = Integer.valueOf(args[1]);

        opticalDups = Integer.valueOf(args[2]);

      }

      catch (NumberFormatException error) {

        System.err.println("When called with three arguments, must be integers");

        System.exit(1);

      }

      readPairs += uniqueReadPairs;

    }

    NumberFormat decimalFormat = NumberFormat.getPercentInstance();

    decimalFormat.setMinimumFractionDigits(2);

    decimalFormat.setMaximumFractionDigits(2);

    long result2 = readPairs*readPairs/(2*(readPairs-uniqueReadPairs));

    System.out.print("Unique Reads: " + NumberFormat.getInstance().format(uniqueReadPairs) + " ");

    if (totalReadPairs > 0) {

      System.out.println("(" + decimalFormat.format(uniqueReadPairs/(double)totalReadPairs) + ")");

    }

    else {

      System.out.println();

    }

    System.out.print("PCR Duplicates: " + nf.format(readPairs - uniqueReadPairs) + " ");

    if (totalReadPairs > 0) {

      System.out.println("(" + decimalFormat.format((readPairs - uniqueReadPairs)/(double)totalReadPairs) + ")");

    }

    else {

      System.out.println();

    }

    System.out.print("Optical Duplicates: " + nf.format(opticalDups) + " ");

    if (totalReadPairs > 0) {

      System.out.println("(" + decimalFormat.format(opticalDups/(double)totalReadPairs) + ")");

    }

    else {

      System.out.println();

    }

    long result;

    try {

      result = estimateLibrarySize(readPairs, uniqueReadPairs);

    }

    catch (NullPointerException e) {

      System.err.println("Library complexity undefined when total = " + readPairs + " and unique = " + uniqueReadPairs);

      return;

    }

    System.out.println("Library Complexity Estimate: " + nf.format(result));

    //    System.out.println("Library complexity (old): " + nf.format(result2));

  }

  /**

   * Estimates the size of a library based on the number of paired end molecules 

   * observed and the number of unique pairs observed.

   * <br>

   * Based on the Lander-Waterman equation that states:<br>

   *     C/X = 1 - exp( -N/X )<br>

   * where<br>

   *     X = number of distinct molecules in library<br>

   *     N = number of read pairs<br>

   *     C = number of distinct fragments observed in read pairs<br>

   */

  public static Long estimateLibrarySize(final long readPairs, 

                                         final long uniqueReadPairs) {

    final long readPairDuplicates = readPairs - uniqueReadPairs;

    if (readPairs > 0 && readPairDuplicates > 0) {

	    long n = readPairs;

	    long c = uniqueReadPairs;

	    double m = 1.0, M = 100.0;

	    if (c >= n || f(m*c, c, n) < 0) {

        throw new IllegalStateException("Invalid values for pairs and unique pairs: " + n + ", " + c);

	    }

	    while( f(M*c, c, n) >= 0 ) {

        m = M;

        M *= 10.0;

	    }

	    double r = (m+M)/2.0;

	    double u = f( r * c, c, n );

	    int i = 0;

	    while (u != 0 && i < 1000) {

        if (u > 0) m = r;

        else M = r;

        r = (m+M)/2.0;

        u = f( r * c, c, n );

        i++;

	    }

	    if (i == 1000) {

        System.err.println("Iterated 1000 times, returning estimate");

	    }

	    return (long) (c * (m+M)/2.0);

    }

    else {

	    return null;

    }

  }

  /** Method that is used in the computation of estimated library size. */

  private static double f(double x, double c, double n) {

    return c/x - 1 + Math.exp(-n/x);

  }

}    

#!/bin/bash

##########

#The MIT License (MIT)

#

# Copyright (c) 2015 Aiden Lab

#

# Permission is hereby granted, free of charge, to any person obtaining a copy

# of this software and associated documentation files (the "Software"), to deal

# in the Software without restriction, including without limitation the rights

# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

# copies of the Software, and to permit persons to whom the Software is

# furnished to do so, subject to the following conditions:

#

# The above copyright notice and this permission notice shall be included in

# all copies or substantial portions of the Software.

#

#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,

#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN

#  THE SOFTWARE.

##########

#

# Sanity check once pipeline is complete to be sure number of lines adds up, deduping

# proceeded appropriately, and that files were created

# Juicer version 1.5

# Start by checking the statistics to see if they add up 

res1=`cat ${splitdir}/*norm*res*`

check1=`cat ${splitdir}/*norm*res* | awk '{s2+=$2; s3+=$3; s4+=$4; s5+=$5; s6+=$6}END{if (s2 != s3+s4+s5+s6){print 0}else{print 1}}'`

if [ $check1 -eq 0 ] || [ -z "$res1" ]

then

    echo "***! Error! The statistics do not add up. Alignment likely failed to complete on one or more files. Run relaunch_prep.sh"

    echo "Stats don't add up.  Check ${outputdir} for results"

    exit 1

fi

# Check the sizes of merged_sort versus the dups/no dups files to be sure

# no reads were lost

total=1

total2=0

total=`ls -l ${outputdir}/merged_sort.txt | awk '{print $5}'`

total2=`ls -l ${outputdir}/merged_nodups.txt ${outputdir}/dups.txt ${outputdir}/opt_dups.txt | awk '{sum = sum + $5}END{print sum}'`

if [ -z $total ] || [ -z $total2 ] || [ $total -ne $total2 ]

then

    echo "***! Error! The sorted file and dups/no dups files do not add up, or were empty. Merge or dedupping likely failed, restart pipeline with -S merge or -S dedup"

    echo "Dups don't add up.  Check ${outputdir} for results"

    exit 1

fi

wctotal=`cat ${splitdir}/*_linecount.txt | awk '{sum+=$1}END{print sum/4}'`

check2=`cat ${splitdir}/*norm*res* | awk '{s2+=$2;}END{print s2}'`

if [ $wctotal -ne $check2 ]

then

    echo "***! Error! The number of reads in the fastqs (${wctotal}) is not the same as the number of reads reported in the stats (${check2}), likely due to a failure during alignment"

    echo "Reads don't add up.  Check ${outputdir} for results"

    exit 1

fi

if [ -n "$early" ]

then

    echo "(-: Pipeline successfully completed (-:";

    echo "Run cleanup.sh to remove the splits directory";

    echo "Check ${outputdir} for results"

elif [ -f ${outputdir}/inter.hic ] && [ -s ${outputdir}/inter.hic ] && [ -f ${outputdir}/inter_30.hic ] && [ -s ${outputdir}/inter_30.hic ]

then

    echo "(-: Pipeline successfully completed (-:";

    echo "Run cleanup.sh to remove the splits directory";

    echo "Check ${outputdir} for results"

else

    echo "***! Error! Either inter.hic or inter_30.hic were not created"

    echo "Either inter.hic or inter_30.hic were not created.  Check ${outputdir} for results"

    exit 1

fi