[med-svn] [gmap] 01/03: Imported Upstream version 2015-11-20
Alex Mestiashvili
malex-guest at moszumanska.debian.org
Thu Nov 26 18:54:54 UTC 2015
This is an automated email from the git hooks/post-receive script.
malex-guest pushed a commit to branch master
in repository gmap.
commit ea9d733e791df53b1ee0f7d815c7ee4f1050fe4d
Author: Alexandre Mestiashvili <alex at biotec.tu-dresden.de>
Date: Thu Nov 26 11:51:04 2015 +0100
Imported Upstream version 2015-11-20
---
ChangeLog | 16 +++
VERSION | 2 +-
configure | 24 ++---
src/bitpack64-write.c | 28 +++---
src/bitpack64-write.h | 14 +--
src/gsnap.c | 6 +-
src/pair.c | 110 ++++++++++++++-------
src/samprint.c | 269 ++++++++++++++++++++++++++++----------------------
src/stage3hr.c | 197 ++++++++++++++++++++++--------------
9 files changed, 404 insertions(+), 262 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index c721c9a..5f82c15 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,19 @@
+2015-11-20 twu
+
+ * VERSION: Updated version number
+
+ * samprint.c: Applied revision 179289 from trunk to fix some issues with
+ --add-paired-nomappers option
+
+ * pair.c, stage3hr.c: Applied revisions 179286 and 179286 from trunk to
+ exclude alignments to circular chromosomes that extend below the first
+ copy or above the second copy
+
+ * gsnap.c: Fixed check on floating point values for --min-coverage
+
+ * bitpack64-write.c, bitpack64-write.h: Changed variable type for n from
+ UINT4 to Oligospace_T, to allow for k-mer of 16
+
2015-09-30 twu
* VERSION: Updated version number
diff --git a/VERSION b/VERSION
index 0189fb2..b8b724d 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-2015-09-29
\ No newline at end of file
+2015-11-20
\ No newline at end of file
diff --git a/configure b/configure
index 9df6267..53d74b8 100755
--- a/configure
+++ b/configure
@@ -1,6 +1,6 @@
#! /bin/sh
# Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.63 for gmap 2015-09-29.
+# Generated by GNU Autoconf 2.63 for gmap 2015-11-20.
#
# Report bugs to <Thomas Wu <twu at gene.com>>.
#
@@ -745,8 +745,8 @@ SHELL=${CONFIG_SHELL-/bin/sh}
# Identity of this package.
PACKAGE_NAME='gmap'
PACKAGE_TARNAME='gmap'
-PACKAGE_VERSION='2015-09-29'
-PACKAGE_STRING='gmap 2015-09-29'
+PACKAGE_VERSION='2015-11-20'
+PACKAGE_STRING='gmap 2015-11-20'
PACKAGE_BUGREPORT='Thomas Wu <twu at gene.com>'
ac_unique_file="src/gmap.c"
@@ -1513,7 +1513,7 @@ if test "$ac_init_help" = "long"; then
# Omit some internal or obsolete options to make the list less imposing.
# This message is too long to be a string in the A/UX 3.1 sh.
cat <<_ACEOF
-\`configure' configures gmap 2015-09-29 to adapt to many kinds of systems.
+\`configure' configures gmap 2015-11-20 to adapt to many kinds of systems.
Usage: $0 [OPTION]... [VAR=VALUE]...
@@ -1584,7 +1584,7 @@ fi
if test -n "$ac_init_help"; then
case $ac_init_help in
- short | recursive ) echo "Configuration of gmap 2015-09-29:";;
+ short | recursive ) echo "Configuration of gmap 2015-11-20:";;
esac
cat <<\_ACEOF
@@ -1721,7 +1721,7 @@ fi
test -n "$ac_init_help" && exit $ac_status
if $ac_init_version; then
cat <<\_ACEOF
-gmap configure 2015-09-29
+gmap configure 2015-11-20
generated by GNU Autoconf 2.63
Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
@@ -1735,7 +1735,7 @@ cat >config.log <<_ACEOF
This file contains any messages produced by compilers while
running configure, to aid debugging if configure makes a mistake.
-It was created by gmap $as_me 2015-09-29, which was
+It was created by gmap $as_me 2015-11-20, which was
generated by GNU Autoconf 2.63. Invocation command line was
$ $0 $@
@@ -2105,8 +2105,8 @@ ac_compiler_gnu=$ac_cv_c_compiler_gnu
{ $as_echo "$as_me:$LINENO: checking package version" >&5
$as_echo_n "checking package version... " >&6; }
-{ $as_echo "$as_me:$LINENO: result: 2015-09-29" >&5
-$as_echo "2015-09-29" >&6; }
+{ $as_echo "$as_me:$LINENO: result: 2015-11-20" >&5
+$as_echo "2015-11-20" >&6; }
### Read defaults
@@ -4172,7 +4172,7 @@ fi
# Define the identity of the package.
PACKAGE='gmap'
- VERSION='2015-09-29'
+ VERSION='2015-11-20'
cat >>confdefs.h <<_ACEOF
@@ -26591,7 +26591,7 @@ exec 6>&1
# report actual input values of CONFIG_FILES etc. instead of their
# values after options handling.
ac_log="
-This file was extended by gmap $as_me 2015-09-29, which was
+This file was extended by gmap $as_me 2015-11-20, which was
generated by GNU Autoconf 2.63. Invocation command line was
CONFIG_FILES = $CONFIG_FILES
@@ -26654,7 +26654,7 @@ Report bugs to <bug-autoconf at gnu.org>."
_ACEOF
cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
ac_cs_version="\\
-gmap config.status 2015-09-29
+gmap config.status 2015-11-20
configured by $0, generated by GNU Autoconf 2.63,
with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\"
diff --git a/src/bitpack64-write.c b/src/bitpack64-write.c
index 6314470..85ca9f8 100644
--- a/src/bitpack64-write.c
+++ b/src/bitpack64-write.c
@@ -1,4 +1,4 @@
-static char rcsid[] = "$Id: bitpack64-write.c 153955 2014-11-24 17:54:45Z twu $";
+static char rcsid[] = "$Id: bitpack64-write.c 179363 2015-11-20 22:13:52Z twu $";
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
@@ -4573,11 +4573,11 @@ compute_q4_diffs_bidir_huge (UINT4 *diffs, UINT8 *values) {
possibly stored as the final metainfo value */
/* Stored in columnar order */
void
-Bitpack64_write_differential (char *ptrsfile, char *compfile, UINT4 *ascending, UINT4 n) {
+Bitpack64_write_differential (char *ptrsfile, char *compfile, UINT4 *ascending, Oligospace_T n) {
FILE *ptrs_fp, *comp_fp;
UINT4 *ptrs;
int ptri, i;
- UINT4 positioni;
+ Oligospace_T positioni;
/* Buffer is used to avoid frequent writes to the file */
UINT4 *buffer;
@@ -4592,6 +4592,8 @@ Bitpack64_write_differential (char *ptrsfile, char *compfile, UINT4 *ascending,
write_setup();
+ printf("Entered Bitpack64_write_differential with n %llu\n",n);
+
/* 2 metavalues: nwritten (pointer) and cumulative sum for block.
Packsize can be computed from difference between successive
pointers, if only even packsizes are allowed */
@@ -4688,11 +4690,11 @@ Bitpack64_write_differential (char *ptrsfile, char *compfile, UINT4 *ascending,
possibly stored as the final metainfo value */
/* D4 stored in columnar order, plus D1 stored as direct */
void
-Bitpack64_write_differential_paired (char *ptrsfile, char *compfile, UINT4 *ascending, UINT4 n) {
+Bitpack64_write_differential_paired (char *ptrsfile, char *compfile, UINT4 *ascending, Oligospace_T n) {
FILE *ptrs_fp, *comp_fp;
UINT4 *ptrs;
int ptri, i;
- UINT4 positioni;
+ Oligospace_T positioni;
/* Buffer is used to avoid frequent writes to the file */
UINT4 *buffer;
@@ -4831,14 +4833,14 @@ Bitpack64_write_differential_paired (char *ptrsfile, char *compfile, UINT4 *asce
possibly stored as the final metainfo value */
/* Stored in columnar order */
void
-Bitpack64_write_fixed10 (char *ptrsfile, char *compfile, UINT4 *ascending, UINT4 n) {
+Bitpack64_write_fixed10 (char *ptrsfile, char *compfile, UINT4 *ascending, Oligospace_T n) {
#ifndef USE_ONE_FILE_FOR_FIXED
FILE *ptrs_fp;
#endif
FILE *comp_fp;
UINT4 *ptrs;
int ptri, i;
- UINT4 positioni;
+ Oligospace_T positioni;
/* Buffer is used to avoid frequent writes to the file */
UINT4 *buffer;
@@ -4969,13 +4971,13 @@ Bitpack64_write_fixed10 (char *ptrsfile, char *compfile, UINT4 *ascending, UINT4
void
Bitpack64_write_differential_huge (char *pagesfile, char *ptrsfile, char *compfile,
- UINT8 *ascending, UINT4 n) {
+ UINT8 *ascending, Oligospace_T n) {
UINT8 currpage, nextpage;
FILE *pages_fp, *ptrs_fp, *comp_fp;
UINT4 pages[25]; /* Allows us to handle up to 100 billion positions */
UINT4 *ptrs;
int ptri;
- UINT4 positioni;
+ Oligospace_T positioni;
/* Buffer is used to avoid frequent writes to the file */
UINT4 *buffer;
@@ -5129,7 +5131,7 @@ Bitpack64_write_differential_huge (char *pagesfile, char *ptrsfile, char *compfi
void
Bitpack64_write_fixed10_huge (char *pagesfile, char *ptrsfile, char *compfile,
- UINT8 *ascending, UINT4 n) {
+ UINT8 *ascending, Oligospace_T n) {
#ifndef USE_ONE_FILE_FOR_FIXED
FILE *ptrs_fp;
#endif
@@ -5138,7 +5140,7 @@ Bitpack64_write_fixed10_huge (char *pagesfile, char *ptrsfile, char *compfile,
UINT4 pages[25]; /* Allows us to handle up to 100 billion positions */
UINT4 *ptrs;
int ptri;
- UINT4 positioni;
+ Oligospace_T positioni;
/* Buffer is used to avoid frequent writes to the file */
UINT4 *buffer;
@@ -5352,11 +5354,11 @@ compute_packsize (UINT4 *values) {
/* Want to store values 0..n-1. The value direct[n] does not exist. */
/* Stored in vertical order */
void
-Bitpack64_write_direct (char *ptrsfile, char *compfile, UINT4 *direct, UINT4 n) {
+Bitpack64_write_direct (char *ptrsfile, char *compfile, UINT4 *direct, Oligospace_T n) {
FILE *ptrs_fp, *comp_fp;
UINT4 *ptrs;
int ptri, i;
- UINT4 positioni;
+ Oligospace_T positioni;
UINT4 *buffer;
int buffer_size = BUFFER_SIZE;
diff --git a/src/bitpack64-write.h b/src/bitpack64-write.h
index 1688c30..693c821 100644
--- a/src/bitpack64-write.h
+++ b/src/bitpack64-write.h
@@ -1,4 +1,4 @@
-/* $Id: bitpack64-write.h 165968 2015-05-20 00:15:38Z twu $ */
+/* $Id: bitpack64-write.h 179363 2015-11-20 22:13:52Z twu $ */
#ifndef BITPACK64_WRITE_INCLUDED
#define BITPACK64_WRITE_INCLUDED
#include <stdio.h>
@@ -6,20 +6,20 @@
/* Stores the $(n+1)$ values [0..n] */
extern void
-Bitpack64_write_differential (char *ptrsfile, char *compfile, UINT4 *ascending, UINT4 n);
+Bitpack64_write_differential (char *ptrsfile, char *compfile, UINT4 *ascending, Oligospace_T n);
extern void
-Bitpack64_write_differential_paired (char *ptrsfile, char *compfile, UINT4 *ascending, UINT4 n);
+Bitpack64_write_differential_paired (char *ptrsfile, char *compfile, UINT4 *ascending, Oligospace_T n);
extern void
-Bitpack64_write_fixed10 (char *ptrsfile, char *compfile, UINT4 *ascending, UINT4 n);
+Bitpack64_write_fixed10 (char *ptrsfile, char *compfile, UINT4 *ascending, Oligospace_T n);
extern void
Bitpack64_write_differential_huge (char *pagesfile, char *ptrsfile, char *compfile,
- UINT8 *ascending, UINT4 n);
+ UINT8 *ascending, Oligospace_T n);
extern void
Bitpack64_write_fixed10_huge (char *pagesfile, char *ptrsfile, char *compfile,
- UINT8 *ascending, UINT4 n);
+ UINT8 *ascending, Oligospace_T n);
/* Stores the $n$ values [0..(n-1)] */
extern void
-Bitpack64_write_direct (char *ptrsfile, char *compfile, UINT4 *direct, UINT4 n);
+Bitpack64_write_direct (char *ptrsfile, char *compfile, UINT4 *direct, Oligospace_T n);
#endif
diff --git a/src/gsnap.c b/src/gsnap.c
index 4660222..e565155 100644
--- a/src/gsnap.c
+++ b/src/gsnap.c
@@ -1,4 +1,4 @@
-static char rcsid[] = "$Id: gsnap.c 173896 2015-09-12 00:11:40Z twu $";
+static char rcsid[] = "$Id: gsnap.c 179364 2015-11-20 22:14:22Z twu $";
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
@@ -1868,10 +1868,6 @@ parse_command_line (int argc, char *argv[], int optind) {
if (user_mincoverage_float > 1.0 && user_mincoverage_float != rint(user_mincoverage_float)) {
fprintf(stderr,"Cannot specify fractional value %f for --max-mismatches except between 0.0 and 1.0\n",user_mincoverage_float);
return 9;
- } else if (user_mincoverage_float > 0.10 && user_mincoverage_float < 1.0) {
- fprintf(stderr,"Your value %f for --max-mismatches implies more than 10%% mismatches, which does not make sense\n",
- user_mincoverage_float);
- return 9;
}
} else if (!strcmp(long_name,"indel-endlength")) {
diff --git a/src/pair.c b/src/pair.c
index 60b6e26..3fd1649 100644
--- a/src/pair.c
+++ b/src/pair.c
@@ -1,4 +1,4 @@
-static char rcsid[] = "$Id: pair.c 174482 2015-09-22 00:58:39Z twu $";
+static char rcsid[] = "$Id: pair.c 179365 2015-11-20 22:15:56Z twu $";
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
@@ -3149,53 +3149,97 @@ Pair_print_gff3 (Filestring_T fp, struct T *pairs, int npairs, int pathnum, char
int
Pair_circularpos (int *alias, struct T *pairs, int npairs, Chrpos_T chrlength, bool plusp, int querylength) {
- int i;
+ Chrpos_T low, high;
struct T *ptr;
+ int i, ninsertions, querypos;
if (plusp == true) {
- i = 0;
- ptr = pairs;
- if (ptr->genomepos >= chrlength) {
+ low = pairs[0].genomepos /*- pairs[0].querypos*/;
+ high = pairs[npairs-1].genomepos /*+ (querylength - pairs[npairs-1].querypos)*/;
+
+ if (low >= chrlength) {
+ /* ptr->genomepos + ptr->querypos >= chrlength (SOFT_CLIPS_AVOID_CIRCULARIZATION */
+ /* ptr->genomepos >= chrlength */
+
/* All of read after trimming is in circular alias */
- *alias = +1;
+ if (high > chrlength + chrlength) {
+ *alias = +2;
+ } else {
+ *alias = +1;
+ }
return -1;
+
+ } else if (high < chrlength) {
+ /* All of read after trimming is in circular proper */
+ if (pairs[0].genomepos < (Chrpos_T) pairs[0].querypos) {
+ *alias = -2;
+ } else {
+ *alias = -1;
+ }
+ return -1;
+
} else {
- while (i < npairs && ptr->genomepos < chrlength) {
- i++;
+ /* Some of read is in circular proper and some is in circular alias */
+ i = 0;
+ ptr = pairs;
+ ninsertions = 0;
+
+ while (i++ < npairs && ptr->genomepos < chrlength) {
+ querypos = ptr->querypos;
+ if (ptr->genome == ' ' && ptr->gapp == false) {
+ ninsertions += 1;
+ }
ptr++;
}
- if (i >= npairs) {
- /* All of read after trimming is in circular proper */
- *alias = -1;
- return -1;
- } else {
- /* Some of read is in circular proper and some is in circular alias */
- *alias = 0;
- return ptr->querypos;
- }
+
+ *alias = 0;
+ return querypos - ninsertions;
}
} else {
- i = npairs - 1;
- ptr = &(pairs[i]);
- if (ptr->genomepos >= chrlength) {
+ low = pairs[npairs-1].genomepos /*- (querylength - pairs[npairs-1].querypos)*/;
+ high = pairs[0].genomepos /*+ pairs[0].querypos*/;
+
+ if (low >= chrlength) {
+ /* ptr->genomepos + (querylength - ptr->querypos) >= chrlength (SOFT_CLIPS_AVOID_CIRCULARIZATION) */
+ /* ptr->genomepos >= chrlength */
+
/* All of read after trimming is in circular alias */
- *alias = +1;
- return -1;
- } else {
- while (i >= 0 && ptr->genomepos < chrlength) {
- i--;
- ptr--;
+ if (high > chrlength + chrlength) {
+ *alias = +2;
+ } else {
+ *alias = +1;
}
- if (i < 0) {
- /* All of read after trimming is in circular proper */
- *alias = -1;
- return -1;
+ return -1;
+
+ } else if (high < chrlength) {
+ /* All of read after trimming is in circular proper */
+ if (pairs[npairs-1].genomepos < (Chrpos_T) (querylength - pairs[npairs-1].querypos)) {
+ *alias = -2;
} else {
- /* Some of read is in circular proper and some is in circular alias */
- *alias = 0;
- return (querylength - ptr->querypos/*- 1*/);
+ *alias = -1;
}
+ return -1;
+
+ } else {
+ /* Some of read is in circular proper and some is in circular alias */
+ i = npairs - 1;
+ ptr = &(pairs[i]);
+ ninsertions = 0;
+
+ while (--i >= 0 && ptr->genomepos < chrlength) {
+ querypos = ptr->querypos;
+ if (ptr->genome == ' ' && ptr->gapp == false) {
+ ninsertions += 1;
+ }
+ --ptr;
+ }
+ /* printf("Ended up at querypos %d, genomepos %u < chrlength %u\n",ptr->querypos,ptr->genomepos,chrlength); */
+ /* printf("Returning querylength %d - querypos %d - ninsertions %d = %d\n",
+ querylength,ptr->querypos,ninsertions,querylength - ptr->querypos - ninsertions); */
+
+ *alias = 0;
+ return (querylength - querypos - ninsertions);
}
}
}
diff --git a/src/samprint.c b/src/samprint.c
index 19a2ca6..27d7c10 100644
--- a/src/samprint.c
+++ b/src/samprint.c
@@ -1,4 +1,4 @@
-static char rcsid[] = "$Id: samprint.c 172734 2015-08-27 16:35:15Z twu $";
+static char rcsid[] = "$Id: samprint.c 179366 2015-11-20 22:16:36Z twu $";
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
@@ -4150,17 +4150,7 @@ SAM_print_paired (Filestring_T fp, Filestring_T fp_failedinput_1, Filestring_T f
resulttype,/*first_read_p*/true,/*npaths_mate*/1,quality_shift,sam_read_group_id,
invert_first_p,invert_second_p,merge_samechr_p);
- if (add_paired_nomappers_p == true) {
- /* matching nomapper for second end */
- SAM_print_nomapping(fp,abbrev,queryseq2,/*mate*/hit5,acc1,acc2,chromosome_iit,
- resulttype,/*first_read_p*/false,/*npaths*/1,/*npaths_mate*/1,/*mate_chrpos*/chrpos5,
- quality_shift,sam_read_group_id,invert_first_p,invert_second_p);
-
- /* matching nomapper for first end */
- SAM_print_nomapping(fp,abbrev,queryseq1,/*mate*/hit3,acc1,acc2,chromosome_iit,
- resulttype,/*first_read_p*/true,/*npaths*/1,/*npaths_mate*/1,/*mate_chrpos*/chrpos3,
- quality_shift,sam_read_group_id,invert_first_p,invert_second_p);
- }
+ /* Note: Do not act on add_paired_nomappers_p, since the two reads are artificially paired up already */
/* print second end */
/* Stage3end_eval_and_sort(stage3array2,npaths2,maxpaths_report,queryseq2); */
@@ -4210,60 +4200,19 @@ SAM_print_paired (Filestring_T fp, Filestring_T fp_failedinput_1, Filestring_T f
}
#endif
- /* print first end results */
- if (npaths2 == 0) {
- mate = (Stage3end_T) NULL;
- chrpos3 = 0U;
- } else if (quiet_if_excessive_p && npaths2 > maxpaths_report) {
- mate = (Stage3end_T) NULL;
- chrpos3 = 0U;
- } else {
- mate = stage3array2[0];
- hardclip3_low = hardclip3_high = 0;
- chrpos3 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,mate,Shortread_fulllength(queryseq2),/*first_read_p*/false);
- }
-
- if (npaths1 == 1) {
- stage3 = stage3array1[0];
- hardclip5_low = hardclip5_high = 0;
- chrpos5 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,stage3,Shortread_fulllength(queryseq1),/*first_read_p*/true);
-
- SAM_print(fp,fp_failedinput_1,abbrev,stage3,mate,acc1,acc2,/*pathnum*/1,npaths1,
- Stage3end_absmq_score(stage3),first_absmq1,second_absmq1,
- Stage3end_mapq_score(stage3),chromosome_iit,
- /*queryseq*/queryseq1,/*queryseq_mate*/queryseq2,
- /*pairedlength*/0U,/*chrpos*/chrpos5,/*mate_chrpos*/chrpos3,
- /*clipdir*/0,/*hardclip5_low*/0,/*hardclip5_high*/0,/*hardclip3_low*/0,/*hardclip3_high*/0,
- resulttype,/*first_read_p*/true,/*npaths_mate*/npaths2,quality_shift,sam_read_group_id,
- invert_first_p,invert_second_p,merge_samechr_p);
+ if (add_paired_nomappers_p == true) {
+ /* Artificially pair up results */
+ for (pathnum = 1; pathnum <= npaths1 && pathnum <= npaths2 && pathnum <= maxpaths_report; pathnum++) {
+ /* hardclip5_low = hardclip5_high = 0; */
+ chrpos5 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,/*stage3*/stage3array1[pathnum-1],
+ Shortread_fulllength(queryseq1),/*first_read_p*/true);
- if (add_paired_nomappers_p == true) {
- /* matching nomapper for second end */
- SAM_print_nomapping(fp,abbrev,queryseq2,/*mate*/stage3,acc1,acc2,chromosome_iit,
- resulttype,/*first_read_p*/false,/*npaths*/1,/*npaths_mate*/npaths2,
- /*mate_chrpos*/chrpos5,
- quality_shift,sam_read_group_id,invert_first_p,invert_second_p);
- }
+ /* hardclip3_low = hardclip3_high = 0; */
+ chrpos3 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,/*stage3*/stage3array2[pathnum-1],
+ Shortread_fulllength(queryseq2),/*first_read_p*/false);
- } else if (quiet_if_excessive_p && npaths1 > maxpaths_report) {
- if (add_paired_nomappers_p == true) {
- /* Handle nomappers with each mapped mate */
- } else {
- /* Just printing one end as nomapping */
- SAM_print_nomapping(fp,abbrev,queryseq1,mate,acc1,acc2,chromosome_iit,
- resulttype,/*first_read_p*/true,npaths1,/*npaths_mate*/npaths2,
- /*mate_chrpos*/chrpos3,
- quality_shift,sam_read_group_id,invert_first_p,invert_second_p);
- }
-
- } else {
- for (pathnum = 1; pathnum <= npaths1 && pathnum <= maxpaths_report; pathnum++) {
stage3 = stage3array1[pathnum-1];
- hardclip5_low = hardclip5_high = 0;
- chrpos5 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/hardclip5_high,stage3,Shortread_fulllength(queryseq1),
- /*first_read_p*/true);
-
- SAM_print(fp,fp_failedinput_1,abbrev,stage3,mate,acc1,acc2,pathnum,npaths1,
+ SAM_print(fp,fp_failedinput_1,abbrev,stage3,/*mate*/stage3array2[pathnum-1],acc1,acc2,pathnum,npaths1,
Stage3end_absmq_score(stage3),first_absmq1,second_absmq1,
Stage3end_mapq_score(stage3),chromosome_iit,
/*queryseq*/queryseq1,/*queryseq_mate*/queryseq2,
@@ -4272,79 +4221,141 @@ SAM_print_paired (Filestring_T fp, Filestring_T fp_failedinput_1, Filestring_T f
resulttype,/*first_read_p*/true,/*npaths_mate*/npaths2,quality_shift,sam_read_group_id,
invert_first_p,invert_second_p,merge_samechr_p);
- if (add_paired_nomappers_p == true) {
+ stage3 = stage3array2[pathnum-1];
+ SAM_print(fp,fp_failedinput_2,abbrev,stage3,/*mate*/stage3array1[pathnum-1],acc1,acc2,pathnum,npaths2,
+ Stage3end_absmq_score(stage3),first_absmq2,second_absmq2,
+ Stage3end_mapq_score(stage3),chromosome_iit,
+ /*queryseq*/queryseq2,/*queryseq_mate*/queryseq1,
+ /*pairedlength*/0U,/*chrpos*/chrpos3,/*mate_chrpos*/chrpos5,
+ /*clipdir*/0,/*hardclip5_low*/0,/*hardclip5_high*/0,/*hardclip3_low*/0,/*hardclip3_high*/0,
+ resulttype,/*first_read_p*/false,/*npaths_mate*/npaths1,quality_shift,sam_read_group_id,
+ invert_second_p,invert_first_p,merge_samechr_p);
+ }
+
+ /* Print remaining results with non-mappers */
+ if (npaths1 > npaths2) {
+ for ( ; pathnum <= npaths1 && pathnum <= maxpaths_report; pathnum++) {
+ stage3 = stage3array1[pathnum-1];
+ /* hardclip5_low = hardclip5_high = 0; */
+ chrpos5 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,stage3,
+ Shortread_fulllength(queryseq1),/*first_read_p*/true);
+ chrpos3 = 0;
+
+ SAM_print(fp,fp_failedinput_1,abbrev,stage3,/*mate*/NULL,acc1,acc2,pathnum,npaths1,
+ Stage3end_absmq_score(stage3),first_absmq1,second_absmq1,
+ Stage3end_mapq_score(stage3),chromosome_iit,
+ /*queryseq*/queryseq1,/*queryseq_mate*/queryseq2,
+ /*pairedlength*/0U,/*chrpos*/chrpos5,/*mate_chrpos*/chrpos3,
+ /*clipdir*/0,/*hardclip5_low*/0,/*hardclip5_high*/0,/*hardclip3_low*/0,/*hardclip3_high*/0,
+ resulttype,/*first_read_p*/true,/*npaths_mate*/npaths2,quality_shift,sam_read_group_id,
+ invert_first_p,invert_second_p,merge_samechr_p);
+
/* matching nomapper for second end */
SAM_print_nomapping(fp,abbrev,queryseq2,/*mate*/stage3,acc1,acc2,chromosome_iit,
resulttype,/*first_read_p*/false,/*npaths*/npaths2,/*npaths_mate*/npaths1,
/*mate_chrpos*/chrpos5,
quality_shift,sam_read_group_id,invert_first_p,invert_second_p);
}
+
+ } else if (npaths2 > npaths1) {
+ for ( ; pathnum <= npaths1 && pathnum <= maxpaths_report; pathnum++) {
+ stage3 = stage3array2[pathnum-1];
+ /* hardclip3_low = hardclip3_high = 0; */
+ chrpos5 = 0;
+ chrpos3 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,stage3,
+ Shortread_fulllength(queryseq2),/*first_read_p*/false);
+
+ /* matching nomapper for first end */
+ SAM_print_nomapping(fp,abbrev,queryseq1,/*mate*/stage3,acc1,acc2,chromosome_iit,
+ resulttype,/*first_read_p*/true,/*npaths*/npaths1,/*npaths_mate*/npaths2,
+ /*mate_chrpos*/chrpos3,
+ quality_shift,sam_read_group_id,invert_first_p,invert_second_p);
+
+ SAM_print(fp,fp_failedinput_2,abbrev,stage3,/*mate*/NULL,acc1,acc2,pathnum,npaths2,
+ Stage3end_absmq_score(stage3),first_absmq2,second_absmq2,
+ Stage3end_mapq_score(stage3),chromosome_iit,
+ /*queryseq*/queryseq2,/*queryseq_mate*/queryseq1,
+ /*pairedlength*/0U,/*chrpos*/chrpos3,/*mate_chrpos*/chrpos5,
+ /*clipdir*/0,/*hardclip5_low*/0,/*hardclip5_high*/0,/*hardclip3_low*/0,/*hardclip3_high*/0,
+ resulttype,/*first_read_p*/false,/*npaths_mate*/npaths1,quality_shift,sam_read_group_id,
+ invert_second_p,invert_first_p,merge_samechr_p);
+ }
}
- }
-
- /* print second end results */
- if (npaths1 == 0) {
- mate = (Stage3end_T) NULL;
- chrpos5 = 0U;
- } else if (quiet_if_excessive_p && npaths1 > maxpaths_report) {
- mate = (Stage3end_T) NULL;
- chrpos5 = 0U;
+
} else {
- mate = stage3array1[0];
- hardclip5_low = hardclip5_high = 0;
- chrpos5 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,mate,Shortread_fulllength(queryseq1),
- /*first_read_p*/true);
- }
+ /* print first end results */
+ if (npaths2 == 0) {
+ mate = (Stage3end_T) NULL;
+ chrpos3 = 0U;
+ } else if (quiet_if_excessive_p && npaths2 > maxpaths_report) {
+ mate = (Stage3end_T) NULL;
+ chrpos3 = 0U;
+ } else {
+ mate = stage3array2[0];
+ hardclip3_low = hardclip3_high = 0;
+ chrpos3 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,mate,Shortread_fulllength(queryseq2),/*first_read_p*/false);
+ }
- if (npaths2 == 1) {
- stage3 = stage3array2[0];
- hardclip3_low = hardclip3_high = 0;
- chrpos3 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,stage3,Shortread_fulllength(queryseq2),
- /*first_read_p*/false);
+ if (npaths1 == 1) {
+ stage3 = stage3array1[0];
+ hardclip5_low = hardclip5_high = 0;
+ chrpos5 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,stage3,Shortread_fulllength(queryseq1),/*first_read_p*/true);
- if (add_paired_nomappers_p == true) {
- /* matching nomapper for first end */
- SAM_print_nomapping(fp,abbrev,queryseq1,/*mate*/stage3,acc1,acc2,chromosome_iit,
- resulttype,/*first_read_p*/true,/*npaths*/npaths1,/*npaths_mate*/npaths2,
+ SAM_print(fp,fp_failedinput_1,abbrev,stage3,mate,acc1,acc2,/*pathnum*/1,npaths1,
+ Stage3end_absmq_score(stage3),first_absmq1,second_absmq1,
+ Stage3end_mapq_score(stage3),chromosome_iit,
+ /*queryseq*/queryseq1,/*queryseq_mate*/queryseq2,
+ /*pairedlength*/0U,/*chrpos*/chrpos5,/*mate_chrpos*/chrpos3,
+ /*clipdir*/0,/*hardclip5_low*/0,/*hardclip5_high*/0,/*hardclip3_low*/0,/*hardclip3_high*/0,
+ resulttype,/*first_read_p*/true,/*npaths_mate*/npaths2,quality_shift,sam_read_group_id,
+ invert_first_p,invert_second_p,merge_samechr_p);
+
+ } else if (quiet_if_excessive_p && npaths1 > maxpaths_report) {
+ /* Just printing one end as nomapping */
+ SAM_print_nomapping(fp,abbrev,queryseq1,mate,acc1,acc2,chromosome_iit,
+ resulttype,/*first_read_p*/true,npaths1,/*npaths_mate*/npaths2,
/*mate_chrpos*/chrpos3,
quality_shift,sam_read_group_id,invert_first_p,invert_second_p);
- }
-
- SAM_print(fp,fp_failedinput_2,abbrev,stage3,mate,acc1,acc2,/*pathnum*/1,npaths2,
- Stage3end_absmq_score(stage3),first_absmq2,second_absmq2,
- Stage3end_mapq_score(stage3),chromosome_iit,
- /*queryseq*/queryseq2,/*queryseq_mate*/queryseq1,
- /*pairedlength*/0U,/*chrpos*/chrpos3,/*mate_chrpos*/chrpos5,
- /*clipdir*/0,/*hardclip5_low*/0,/*hardclip5_high*/0,/*hardclip3_low*/0,/*hardclip3_high*/0,
- resulttype,/*first_read_p*/false,/*npaths_mate*/npaths1,quality_shift,sam_read_group_id,
- invert_second_p,invert_first_p,merge_samechr_p);
- } else if (quiet_if_excessive_p && npaths2 > maxpaths_report) {
- if (add_paired_nomappers_p == true) {
- /* Handle nomappers with each mapped mate */
} else {
- /* Just printing one end as nomapping */
- SAM_print_nomapping(fp,abbrev,queryseq2,mate,acc1,acc2,chromosome_iit,
- resulttype,/*first_read_p*/false,npaths2,/*npaths_mate*/npaths1,
- /*mate_chrpos*/chrpos5,
- quality_shift,sam_read_group_id,invert_second_p,invert_first_p);
+ for (pathnum = 1; pathnum <= npaths1 && pathnum <= maxpaths_report; pathnum++) {
+ stage3 = stage3array1[pathnum-1];
+ hardclip5_low = hardclip5_high = 0;
+ chrpos5 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/hardclip5_high,stage3,Shortread_fulllength(queryseq1),
+ /*first_read_p*/true);
+
+ SAM_print(fp,fp_failedinput_1,abbrev,stage3,mate,acc1,acc2,pathnum,npaths1,
+ Stage3end_absmq_score(stage3),first_absmq1,second_absmq1,
+ Stage3end_mapq_score(stage3),chromosome_iit,
+ /*queryseq*/queryseq1,/*queryseq_mate*/queryseq2,
+ /*pairedlength*/0U,/*chrpos*/chrpos5,/*mate_chrpos*/chrpos3,
+ /*clipdir*/0,/*hardclip5_low*/0,/*hardclip5_high*/0,/*hardclip3_low*/0,/*hardclip3_high*/0,
+ resulttype,/*first_read_p*/true,/*npaths_mate*/npaths2,quality_shift,sam_read_group_id,
+ invert_first_p,invert_second_p,merge_samechr_p);
+ }
+ }
+
+ /* print second end results */
+ if (npaths1 == 0) {
+ mate = (Stage3end_T) NULL;
+ chrpos5 = 0U;
+ } else if (quiet_if_excessive_p && npaths1 > maxpaths_report) {
+ mate = (Stage3end_T) NULL;
+ chrpos5 = 0U;
+ } else {
+ mate = stage3array1[0];
+ hardclip5_low = hardclip5_high = 0;
+ chrpos5 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,mate,Shortread_fulllength(queryseq1),
+ /*first_read_p*/true);
}
- } else {
- for (pathnum = 1; pathnum <= npaths2 && pathnum <= maxpaths_report; pathnum++) {
- stage3 = stage3array2[pathnum-1];
+ if (npaths2 == 1) {
+ stage3 = stage3array2[0];
hardclip3_low = hardclip3_high = 0;
chrpos3 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,stage3,Shortread_fulllength(queryseq2),
/*first_read_p*/false);
-
- if (add_paired_nomappers_p == true) {
- /* matching nomapper for first end */
- SAM_print_nomapping(fp,abbrev,queryseq1,/*mate*/stage3,acc1,acc2,chromosome_iit,
- resulttype,/*first_read_p*/true,/*npaths*/npaths1,/*npaths_mate*/npaths2,
- /*mate_chrpos*/chrpos3,
- quality_shift,sam_read_group_id,invert_first_p,invert_second_p);
- }
- SAM_print(fp,fp_failedinput_2,abbrev,stage3,mate,acc1,acc2,pathnum,npaths2,
+
+ SAM_print(fp,fp_failedinput_2,abbrev,stage3,mate,acc1,acc2,/*pathnum*/1,npaths2,
Stage3end_absmq_score(stage3),first_absmq2,second_absmq2,
Stage3end_mapq_score(stage3),chromosome_iit,
/*queryseq*/queryseq2,/*queryseq_mate*/queryseq1,
@@ -4352,6 +4363,30 @@ SAM_print_paired (Filestring_T fp, Filestring_T fp_failedinput_1, Filestring_T f
/*clipdir*/0,/*hardclip5_low*/0,/*hardclip5_high*/0,/*hardclip3_low*/0,/*hardclip3_high*/0,
resulttype,/*first_read_p*/false,/*npaths_mate*/npaths1,quality_shift,sam_read_group_id,
invert_second_p,invert_first_p,merge_samechr_p);
+
+ } else if (quiet_if_excessive_p && npaths2 > maxpaths_report) {
+ /* Just printing one end as nomapping */
+ SAM_print_nomapping(fp,abbrev,queryseq2,mate,acc1,acc2,chromosome_iit,
+ resulttype,/*first_read_p*/false,npaths2,/*npaths_mate*/npaths1,
+ /*mate_chrpos*/chrpos5,
+ quality_shift,sam_read_group_id,invert_second_p,invert_first_p);
+
+ } else {
+ for (pathnum = 1; pathnum <= npaths2 && pathnum <= maxpaths_report; pathnum++) {
+ stage3 = stage3array2[pathnum-1];
+ hardclip3_low = hardclip3_high = 0;
+ chrpos3 = SAM_compute_chrpos(/*hardclip_low*/0,/*hardclip_high*/0,stage3,Shortread_fulllength(queryseq2),
+ /*first_read_p*/false);
+
+ SAM_print(fp,fp_failedinput_2,abbrev,stage3,mate,acc1,acc2,pathnum,npaths2,
+ Stage3end_absmq_score(stage3),first_absmq2,second_absmq2,
+ Stage3end_mapq_score(stage3),chromosome_iit,
+ /*queryseq*/queryseq2,/*queryseq_mate*/queryseq1,
+ /*pairedlength*/0U,/*chrpos*/chrpos3,/*mate_chrpos*/chrpos5,
+ /*clipdir*/0,/*hardclip5_low*/0,/*hardclip5_high*/0,/*hardclip3_low*/0,/*hardclip3_high*/0,
+ resulttype,/*first_read_p*/false,/*npaths_mate*/npaths1,quality_shift,sam_read_group_id,
+ invert_second_p,invert_first_p,merge_samechr_p);
+ }
}
}
diff --git a/src/stage3hr.c b/src/stage3hr.c
index c62fb51..15a83a3 100644
--- a/src/stage3hr.c
+++ b/src/stage3hr.c
@@ -1,4 +1,4 @@
-static char rcsid[] = "$Id: stage3hr.c 174482 2015-09-22 00:58:39Z twu $";
+static char rcsid[] = "$Id: stage3hr.c 179365 2015-11-20 22:15:56Z twu $";
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
@@ -38,8 +38,9 @@ static char rcsid[] = "$Id: stage3hr.c 174482 2015-09-22 00:58:39Z twu $";
/* Originally added to avoid CIGAR strings like 1S99H, but results in
errors if first chromosome is circular. Now checking in samprint.c
whether the CIGAR string is bad. But now needed again because we
- are allowing alignments that are out of bounds. */
-#define SOFT_CLIPS_AVOID_CIRCULARIZATION 1
+ are allowing alignments that are out of bounds. However, can cause
+ left to be negative, so cannot use. */
+/* #define SOFT_CLIPS_AVOID_CIRCULARIZATION 1 -- Turn off here, but on in pair.c */
#define TRANSLOC_SPECIAL 1 /* Eliminates translocations if non-translocations are found */
@@ -463,7 +464,8 @@ struct T {
bool paired_seenp; /* for paired-end. set to true by Stage3_pair_up(). */
bool concordantp; /* for paired-end. set to true by Stage3_pair_up(). */
- int alias; /* -1 if below chrlength, 0 if straddles or NA (e.g., transloc), and +1 if above */
+ int alias; /* -1 if all below chrlength, 0 if straddles or NA (e.g., transloc), and +1 if above */
+ /* -2 if extends below beginning of circular chromosome, +2 if extends beyond end of second copy */
int circularpos; /* if alias == 0, then amount of queryseq below chrlength */
};
@@ -4580,28 +4582,30 @@ compute_circularpos (int *alias, T hit) {
hit->plusp,hit->querylength);
} else if (hit->plusp == true) {
- if (
-#ifdef SOFT_CLIPS_AVOID_CIRCULARIZATION
- hit->low + hit->trim_left >= hit->chroffset + hit->chrlength
-#else
- hit->low >= hit->chroffset + hit->chrlength
-#endif
- ) {
+ if (hit->low >= hit->chroffset + hit->chrlength) {
+ /* hit->low + hit->trim_left >= hit->chroffset + hit->chrlength (SOFT_CLIPS_AVOID_CIRCULARIZATION) */
+ /* hit->low >= hit->chroffset + hit->chrlength */
+
/* All of read after trimming is in circular alias */
debug12(printf("Soft clip of %d on left avoids circularization\n",hit->trim_left));
- *alias = +1;
+ if (hit->high >= hit->chrhigh) {
+ *alias = +2;
+ } else {
+ *alias = +1;
+ }
return -1;
- } else if (
-#ifdef SOFT_CLIPS_AVOID_CIRCULARIZATION
- hit->high - hit->trim_right <= hit->chroffset + hit->chrlength
-#else
- hit->high <= hit->chroffset + hit->chrlength
-#endif
- ) {
+ } else if (hit->high < hit->chroffset + hit->chrlength) {
+ /* hit->high - hit->trim_right <= hit->chroffset + hit->chrlength (SOFT_CLIPS_AVOID_CIRCULARIZATION) */
+ /* hit->high <= hit->chroffset + hit->chrlength */
+
/* All of read after trimming is in circular proper */
debug12(printf("Soft clip of %d on right avoids circularization\n",hit->trim_right));
- *alias = -1;
+ if (hit->low < hit->chroffset) {
+ *alias = -2;
+ } else {
+ *alias = -1;
+ }
return -1;
} else {
@@ -4617,30 +4621,33 @@ compute_circularpos (int *alias, T hit) {
}
} else {
- if (
-#ifdef SOFT_CLIPS_AVOID_CIRCULARIZATION
- hit->low + hit->trim_right >= hit->chroffset + hit->chrlength
-#else
- hit->low >= hit->chroffset + hit->chrlength
-#endif
- ) {
+ /* Without SOFT_CLIPS_AVOID_CIRCULARIZATION, the branches are similar for plus/minus hits */
+ if (hit->low >= hit->chroffset + hit->chrlength) {
+ /* hit->low + hit->trim_right >= hit->chroffset + hit->chrlength (SOFT_CLIPS_AVOID_CIRCULARIZATION) */
+ /* hit->low >= hit->chroffset + hit->chrlength */
+
/* All of read after trimming is in circular alias */
debug12(printf("Soft clip of %d on right avoids circularization\n",hit->trim_right));
debug12(printf("All of read after trimming is in circular alias\n"));
- *alias = +1;
+ if (hit->high >= hit->chrhigh) {
+ *alias = +2;
+ } else {
+ *alias = +1;
+ }
return -1;
- } else if (
-#ifdef SOFT_CLIPS_AVOID_CIRCULARIZATION
- hit->high - hit->trim_left <= hit->chroffset + hit->chrlength
-#else
- hit->high <= hit->chroffset + hit->chrlength
-#endif
- ) {
+ } else if (hit->high < hit->chroffset + hit->chrlength) {
+ /* hit->high - hit->trim_left <= hit->chroffset + hit->chrlength (SOFT_CLIPS_AVOID_CIRCULARIZATION) */
+ /* hit->high <= hit->chroffset + hit->chrlength */
+
/* All of read after trimming is in circular proper */
debug12(printf("Soft clip of %d on left avoids circularization\n",hit->trim_left));
debug12(printf("All of read after trimming is in circular proper\n"));
- *alias = -1;
+ if (hit->low < hit->chroffset) {
+ *alias = -2;
+ } else {
+ *alias = -1;
+ }
return -1;
} else {
@@ -5131,9 +5138,14 @@ Stage3end_new_substrings (int *found_score, Intlist_T endpoints,
new->concordantp = false;
new->circularpos = compute_circularpos(&new->alias,new);
-
- debug0(printf("Returning %p from Stage3end_new_substrings with found_score %d\n",new,*found_score));
- return new;
+ if (new->alias == +2 || new->alias == -2) {
+ debug0(printf("Returning NULL from Stage3end_new_substrings because of alias of %d\n",new->alias));
+ Stage3end_free(&new);
+ return (T) NULL;
+ } else {
+ debug0(printf("Returning %p from Stage3end_new_substrings with found_score %d\n",new,*found_score));
+ return new;
+ }
}
@@ -5814,8 +5826,12 @@ Stage3end_new_exact (int *found_score, Univcoord_T left, int genomiclength, Comp
new->concordantp = false;
new->circularpos = compute_circularpos(&new->alias,new);
-
- return new;
+ if (new->alias == +2 || new->alias == -2) {
+ Stage3end_free(&new);
+ return (T) NULL;
+ } else {
+ return new;
+ }
}
}
@@ -6014,9 +6030,13 @@ Stage3end_new_substitution (int *found_score, int nmismatches_whole, Univcoord_T
new->concordantp = false;
new->circularpos = compute_circularpos(&new->alias,new);
-
- debug(printf("Returning substitution %p\n",new));
- return new;
+ if (new->alias == +2 || new->alias == -2) {
+ Stage3end_free(&new);
+ return (T) NULL;
+ } else {
+ debug(printf("Returning substitution %p\n",new));
+ return new;
+ }
}
}
@@ -6302,8 +6322,12 @@ Stage3end_new_insertion (int *found_score, int nindels, int indel_pos, int nmism
new->concordantp = false;
new->circularpos = compute_circularpos(&new->alias,new);
-
- return new;
+ if (new->alias == +2 || new->alias == -2) {
+ Stage3end_free(&new);
+ return (T) NULL;
+ } else {
+ return new;
+ }
}
}
@@ -6409,7 +6433,12 @@ Stage3end_new_deletion (int *found_score, int nindels, int indel_pos, int nmisma
alignstart1 = left + (querylength /*- querystart (0)*/);
alignend1 = left + (querylength - indel_pos);
alignstart2 = (left - nindels) + (querylength - indel_pos);
- alignend2 = (left - nindels) /*+ querylength - queryend (querylength)*/;
+ if (left < nindels) {
+ /* Handles deletion at start of genome, and avoids left alignend2 becoming negative */
+ alignend2 = 0;
+ } else {
+ alignend2 = (left - nindels) /*+ querylength - queryend (querylength)*/;
+ }
genomicstart1 = alignstart1;
genomicend1 = alignend1;
@@ -6609,8 +6638,12 @@ Stage3end_new_deletion (int *found_score, int nindels, int indel_pos, int nmisma
new->concordantp = false;
new->circularpos = compute_circularpos(&new->alias,new);
-
- return new;
+ if (new->alias == +2 || new->alias == -2) {
+ Stage3end_free(&new);
+ return (T) NULL;
+ } else {
+ return new;
+ }
}
}
@@ -7078,16 +7111,19 @@ Stage3end_new_splice (int *found_score, int nmismatches_donor, int nmismatches_a
new->concordantp = false;
new->circularpos = compute_circularpos(&new->alias,new);
-
- debug0(printf("Returning new splice %p at genomic %u..%u, donor %p (%u => %u), acceptor %p (%u => %u)\n",
- new,new->genomicstart - new->chroffset,new->genomicend - new->chroffset,donor,
- donor == NULL ? 0 : Substring_left_genomicseg(donor),
- donor == NULL ? 0 : Substring_splicecoord(donor),
- acceptor,acceptor == NULL ? 0 : Substring_left_genomicseg(acceptor),
- acceptor == NULL ? 0 : Substring_splicecoord(acceptor)));
-
- debug0(printf("sensedir %d\n",new->sensedir));
- return new;
+ if (new->alias == +2 || new->alias == -2) {
+ Stage3end_free(&new);
+ return (T) NULL;
+ } else {
+ debug0(printf("Returning new splice %p at genomic %u..%u, donor %p (%u => %u), acceptor %p (%u => %u)\n",
+ new,new->genomicstart - new->chroffset,new->genomicend - new->chroffset,donor,
+ donor == NULL ? 0 : Substring_left_genomicseg(donor),
+ donor == NULL ? 0 : Substring_splicecoord(donor),
+ acceptor,acceptor == NULL ? 0 : Substring_left_genomicseg(acceptor),
+ acceptor == NULL ? 0 : Substring_splicecoord(acceptor)));
+ debug0(printf("sensedir %d\n",new->sensedir));
+ return new;
+ }
}
@@ -7414,8 +7450,12 @@ Stage3end_new_shortexon (int *found_score, Substring_T donor, Substring_T accept
new->concordantp = false;
new->circularpos = compute_circularpos(&new->alias,new);
-
- return new;
+ if (new->alias == +2 || new->alias == -2) {
+ Stage3end_free(&new);
+ return (T) NULL;
+ } else {
+ return new;
+ }
}
@@ -7690,8 +7730,12 @@ Stage3end_new_terminal (int querystart, int queryend, Univcoord_T left, Compress
new->concordantp = false;
new->circularpos = compute_circularpos(&new->alias,new);
-
- return new;
+ if (new->alias == +2 || new->alias == -2) {
+ Stage3end_free(&new);
+ return (T) NULL;
+ } else {
+ return new;
+ }
}
@@ -7996,8 +8040,15 @@ Stage3end_new_gmap (int nmismatches_whole, int nmatches_posttrim, int max_match_
new->concordantp = false;
new->circularpos = compute_circularpos(&new->alias,new);
-
- return new;
+ if (new->alias == +2 || new->alias == -2) {
+ /* Stage3end_free(&new); -- Cannot use, because it frees pairarray */
+ Pair_tokens_free(&new->cigar_tokens);
+ /* No substrings or junctions */
+ FREE(new);
+ return (T) NULL;
+ } else {
+ return new;
+ }
}
@@ -9213,7 +9264,7 @@ Stage3end_noptimal (List_T hitlist, int querylength) {
static Univcoord_T
normalize_coord (Univcoord_T orig, int alias, Chrpos_T chrlength) {
- if (alias > 0) {
+ if (alias == +1) {
return orig - chrlength;
} else {
return orig;
@@ -13147,7 +13198,7 @@ Stage3pair_new (T hit5, T hit3, Univcoord_T *splicesites,
new->insertlength -= hit5->chrlength;
}
- if (hit5->alias > 0) {
+ if (hit5->alias == +1) {
debug0(printf("Unaliasing 5' end\n"));
if (private5p == false) {
new->hit5 = Stage3end_copy(hit5);
@@ -13156,7 +13207,7 @@ Stage3pair_new (T hit5, T hit3, Univcoord_T *splicesites,
unalias_circular(new->hit5);
}
- if (hit3->alias > 0) {
+ if (hit3->alias == +1) {
debug0(printf("Unaliasing 3' end\n"));
if (private3p == false) {
new->hit3 = Stage3end_copy(hit3);
@@ -15408,13 +15459,12 @@ Stage3end_linearize_5 (List_T hitlist) {
if (hit->alias == 0) {
/* Skip */
- } else if (hit->alias > 0) {
+ } else if (hit->alias == +1) {
if (hit->plusp == true) {
unalias_circular(hit);
}
- } else {
- /* hit->alias < 0 */
+ } else if (hit->alias == -1) {
if (hit->plusp == false) {
alias_circular(hit);
}
@@ -15441,13 +15491,12 @@ Stage3end_linearize_3 (List_T hitlist) {
if (hit->alias == 0) {
/* Skip */
- } else if (hit->alias < 0) {
+ } else if (hit->alias == -1) {
if (hit->plusp == true) {
alias_circular(hit);
}
- } else {
- /* hit->alias > 0 */
+ } else if (hit->alias == +1) {
if (hit->plusp == false) {
unalias_circular(hit);
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/gmap.git
More information about the debian-med-commit
mailing list