rollback

#!/usr/bin/awk -f

##########

#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.

##########

# Takes a SAM file, looks for interesting abnormal chimeric reads.

# only those reads mapped to chromosomes 1-24 with MAPQ >= 10, and

# with at least a triple where the minimum distance >= 20Kb

# Juicer version 1.5

function abs(v) {

    return v<0?-v:v;

}

BEGIN{

  OFS="\t";

  mapq_thresh=-1;

  tottot = -1; # will count first non-group

}

{

  # input file is sorted by read name.  Look at read name to group 

  # appropriately

  split($1,a,"/");

  if(a[1]==prev){

    # move on to next record.  look below this block for actions that occur

    # regardless.

    count++;

  }

  else {

    # deal with read pair group

    tottot++;

    printme = 1;

    j = 0;

    for (i in c) {

      split(c[i], tmp);

      mapped[j] = tmp[3] ~ /^[1-9,X,Y,M][0-9,T]?$/ || tmp[3] ~ /^chr[1-9,X,Y,M][0-9,T]?$/;

      name[j] = tmp[1];

      str[j] = tmp[2];

      chr[j] = tmp[3];

      pos[j] = tmp[4];

      m[j] = tmp[5];

      cig[j] = tmp[6];

      seq[j] = tmp[10];

      j++;

    }

    len = j;

    bigdist = 0;

    for (j=0; j<len; j++) {

        printme = printme && mapped[j] && (m[j] >= 10);

        for (k=j+1; k<len; k++) {

            if (abs(pos[j]-pos[k]) >= 20000) {

                bigdist++;

            }

        }

    }

#    printme = printme && (bigdist >= 3);

    if (printme) {

      for (j = 0; j < len; j++) {

          printf("%s %d %s %d %d %s %s ", name[j], str[j], chr[j], pos[j], m[j], cig[j], seq[j]);

      }

      printf("\n");

    }

    # reset variables

    delete c;

    count=1;

    prev=a[1];

  }

  # these happen no matter what, after the above processing

  c[count] = $0;

}

END {

    printme = 1;

    j = 0;

    for (i in c) {

      split(c[i], tmp);

      mapped[j] = tmp[3] ~ /^[1-9,X,Y,M][0-9,T]?$/ || tmp[3] ~ /^chr[1-9,X,Y,M][0-9,T]?$/;

      name[j] = tmp[1];

      str[j] = tmp[2];

      chr[j] = tmp[3];

      pos[j] = tmp[4];

      m[j] = tmp[5];

      cig[j] = tmp[6];

      seq[j] = tmp[10];

      j++;

    }

    len = j;

    bigdist = 0;

    for (j=0; j<len; j++) {

        printme = printme && mapped[j] && (m[j] >= 10);

        for (k=j+1; k<len; k++) {

            if (abs(pos[j]-pos[k]) >= 20000) {

                bigdist++;

            }

        }

    }

#    printme = printme && (bigdist >= 3);

    if (printme) {

      for (j = 0; j < len; j++) {

          printf("%s %d %s %d %d %s %s ", name[j], str[j], chr[j], pos[j], m[j], cig[j], seq[j]);

      }

      printf("\n");

    }

}