[med-svn] [subread] 03/14: Imported Upstream version 1.5.0-p2+dfsg
Alex Mestiashvili
malex-guest at moszumanska.debian.org
Sat Apr 23 09:38:35 UTC 2016
This is an automated email from the git hooks/post-receive script.
malex-guest pushed a commit to branch master
in repository subread.
commit 8b7078e38a46c17b7eb10262470a5a23b1fd0e58
Author: Alexandre Mestiashvili <alex at biotec.tu-dresden.de>
Date: Fri Apr 22 15:53:25 2016 +0200
Imported Upstream version 1.5.0-p2+dfsg
---
doc/SubreadUsersGuide.tex | 20 +-
src/HelperFunctions.c | 148 +++
src/HelperFunctions.h | 1 +
src/Makefile | 2 +-
src/Makefile.Linux | 8 +-
src/SNPCalling.c | 20 +-
src/core-bigtable.c | 40 +-
src/core-indel.c | 19 +-
src/core-interface-aligner.c | 2 +
src/core-interface-subjunc.c | 20 +-
src/core-junction.c | 1166 ++++++++++++++++----
src/core.c | 222 ++--
src/core.h | 3 +-
src/coverage_calc.c | 4 +-
src/gene-algorithms.c | 20 +-
src/gene-algorithms.h | 2 +-
src/gene-value-index.c | 16 +-
src/global-reassembly.c | 5 +-
src/index-builder.c | 2 +
src/input-files.c | 62 +-
src/input-files.h | 5 +-
src/makefile.version | 6 +-
src/propmapped.c | 8 +-
src/read-repair.c | 6 +-
src/readSummary.c | 282 +++--
src/removeDupReads.c | 6 +-
src/sambam-file.c | 6 +
src/sambam-file.h | 1 +
src/subread.h | 2 +-
src/subtools.c | 5 +-
test/featureCounts/data/corner-EXON-ONLY.ora | 9 +
.../data/test-minimum-dup.ora.summary | 12 -
test/featureCounts/test_corner_cases.sh | 1 +
33 files changed, 1668 insertions(+), 463 deletions(-)
diff --git a/doc/SubreadUsersGuide.tex b/doc/SubreadUsersGuide.tex
index c95eee8..48ee9a6 100644
--- a/doc/SubreadUsersGuide.tex
+++ b/doc/SubreadUsersGuide.tex
@@ -35,9 +35,9 @@
\begin{center}
{\Huge\bf Subread/Rsubread Users Guide}\\
\vspace{1 cm}
-{\centering\large Subread v1.5.0/Rsubread v1.20.2\\}
+{\centering\large Subread v1.5.0-p2/Rsubread v1.20.5\\}
\vspace{1 cm}
-\centering 16 December 2015\\
+\centering 11 April 2016\\
\vspace{5 cm}
\Large Wei Shi and Yang Liao\\
\vspace{1 cm}
@@ -47,7 +47,7 @@ The Walter and Eliza Hall Institute of Medical Research\\
The University of Melbourne\\
Melbourne, Australia\\}
\vspace{7 cm}
-\centering Copyright \small{\copyright} 2011 - 2015\\
+\centering Copyright \small{\copyright} 2011 - 2016\\
\end{center}
\end{titlepage}
@@ -419,7 +419,7 @@ Arguments & Description \\
\hline
-b \newline (\code{color2base=TRUE}) & Output base-space reads instead of color-space reads in mapping output for color space data (eg. LifTech SOLiD data). Note that the mapping itself will still be performed at color-space.\\
\hline
--B $<int>$ \newline (\code{nBestLocations}) & Specify the maximal number of equally-best mapping locations allowed to be reported for each read. 1 by default. Allowed values are between 1 to 16 (inclusive). `NH' tag is used to indicate how many alignments are reported for the read and `HI' tag is used for numbering the alignments reported for the same read, in the output. Note that \code{-u} option takes precedence over \code{-B}.\\
+-B $<int>$ \newline (\code{nBestLocations}) & Specify the maximal number of equally-best mapping locations allowed to be reported for each read. 1 by default. `NH' tag is used to indicate how many alignments are reported for the read and `HI' tag is used for numbering the alignments reported for the same read, in the output. Note that \code{-u} option takes precedence over \code{-B}.\\
\hline
-d $<int>$ \newline (\code{minFragLength}) & Specify the minimum fragment/template length, 50 by default. Note that if the two reads from the same pair do not satisfy the fragment length criteria, they will be mapped individually as if they were single-end reads.\\
\hline
@@ -447,7 +447,7 @@ Arguments & Description \\
\hline
-S $<ff:fr:rf>$ \newline (\code{PE\_orientation}) & Specify the orientation of the two reads from the same pair. It has three possible values including `fr', `ff' and `'rf. Letter `f' denotes the forward strand and letter `r' the reverse strand. `fr' by default (ie. the first read in the pair is on the forward strand and the second read on the reverse strand).\\
\hline
-$^*$ -t $<int>$ \newline (\code{type}) & Specify the type of input sequencing data. Possible values include \code{0}, denoting RNA-seq data, or \code{1}, denoting genomic DNA-seq data. For genomic DNA-seq data, the aligner takes into account both the number of matched bases and the number of mis-matched bases to determine the the best mapping location after applying the `seed-and-vote' approach for read mapping. For RNA-seq data, only the number of mis-matched bases is considered for det [...]
+$^*$ -t $<int>$ \newline (\code{type}) & Specify the type of input sequencing data. Possible values include \code{0}, denoting RNA-seq data, or \code{1}, denoting genomic DNA-seq data. Character values including `rna' and `dna' can also be used in the {\R} function. For genomic DNA-seq data, the aligner takes into account both the number of matched bases and the number of mis-matched bases to determine the the best mapping location after applying the `seed-and-vote' approach for read ma [...]
\hline
-T $<int>$ \newline (\code{nthreads}) & Specify the number of threads/CPUs used for mapping. The value should be between 1 and 32. 1 by default.\\
\hline
@@ -841,7 +841,7 @@ input\_files \newline (\code{files}) & Give the names of input read files that i
\hline
-M \newline (\code{countMultiMappingReads}) & If specified, multi-mapping reads/fragments will be counted. A multi-mapping read will be counted up to N times if it has N reported mapping locations. The program uses the `NH' tag to find multi-mapping reads.\\
\hline
--o $<input>$ & Give the name of the output file. The output file contains the number of reads assigned to each meta-feature (or each feature if \code{-f} is specified). Note that the {\featureCounts} function in {\Rsubread} does not use this parameter. It returns a \code{list} object including read summarization results and other data. \\
+-o $<output>$ & Give the name of the output file. The output file contains the number of reads assigned to each meta-feature (or each feature if \code{-f} is specified). Note that the {\featureCounts} function in {\Rsubread} does not use this parameter. It returns a \code{list} object including read summarization results and other data. \\
\hline
-O \newline (\code{allowMultiOverlap}) & If specified, reads (or fragments if \code{-p} is specified) will be allowed to be assigned to more than one matched meta-feature (or feature if \code{-f} is specified). Reads/fragments overlapping with more than one meta-feature/feature will be counted more than once. Note that when performing meta-feature level summarization, a read (or fragment) will still be counted once if it overlaps with multiple features belonging to the same meta-feature [...]
\hline
@@ -855,15 +855,15 @@ input\_files \newline (\code{files}) & Give the names of input read files that i
\hline
-s $<int>$ \newline (\code{isStrandSpecific}) & Indicate if strand-specific read counting should be performed. It has three possible values: 0 (unstranded), 1 (stranded) and 2 (reversely stranded). 0 by default. For paired-end reads, strand of the first read is taken as the strand of the whole fragment and FLAG field of the current read is used to tell if it is the first read in the fragment.\\
\hline
--S $<ff:fr:rf>$ \newline (\code{PE\_orientation}) & Specify the orientation of the two reads from the same pair. It has three possible values including `fr', `ff' and `'rf. Letter `f' denotes the forward strand and letter `r' the reverse strand. `fr' by default (ie. the first read in the pair is on the forward strand and the second read on the reverse strand).\\
-\hline
-t $<input>$ \newline (\code{GTF.featureType}) & Specify the feature type. Only rows which have the matched feature type in the provided GTF annotation file will be included for read counting. `exon' by default.\\
\hline
-T $<int>$ \newline (\code{nthreads}) & Number of the threads. The value should be between 1 and 32. 1 by default.\\
\hline
-v & Output version of the program. \\
\hline
-$--$countSplitAlignmentsOnly \newline (\code{countSplitAlignmentsOnly}) & If specified, only split alignments (CIGAR strings containing letter `N') will be counted. All the other alignments will be ignored. An example of split alignments is the exon-spanning reads in RNA-seq data. If exon-spanning reads need to be assigned to all their overlapping exons, `-f' and `-O' options should be provided as well.\\
+$--$countSplitAlignmentsOnly \newline (\code{splitOnly}) & If specified, only split alignments (CIGAR strings contain letter `N') will be counted. All the other alignments will be ignored. An example of split alignments is the exon-spanning reads in RNA-seq data. If exon-spanning reads need to be assigned to all their overlapping exons, `-f' and `-O' options should be provided as well.\\
+\hline
+$--$countNonSplitAlignmentsOnly \newline (\code{nonSplitOnly}) & If specified, only non-split alignments (CIGAR strings do not contain letter `N') will be counted. All the other alignments will be ignored.\\
\hline
$--$donotsort \newline (\code{autosort}) & If specified, paired end reads will not be re-ordered even if reads from the same pair were found not to be next to each other in the input.\\
\hline
@@ -877,7 +877,7 @@ $--$maxMOp $<int>$ \newline (\code{maxMOp}) & Specify the maximum number of `M'
\hline
$--$minOverlap $<int>$ \newline (\code{minOverlap}) & Specify the minimum required number of overlapping bases between a read (or a fragment) and an overlapping feature. 1 by default. If a negative value is provided, the read will be extended from both ends.\\
\hline
-$--$primary \newline (\code{countPrimaryAlignmentsOnly}) & If specified, only primary alignments will be counted. Primary and secondary alignments are identified using bit 0x100 in the Flag field of SAM/BAM files. All primary alignments in a dataset will be counted no matter they are from multi-mapping reads or not (ie. `-M' is ignored).\\
+$--$primary \newline (\code{primaryOnly}) & If specified, only primary alignments will be counted. Primary and secondary alignments are identified using bit 0x100 in the Flag field of SAM/BAM files. All primary alignments in a dataset will be counted no matter they are from multi-mapping reads or not (ie. `-M' is ignored).\\
\hline
$--$read2pos $<int>$ \newline (\code{read2pos}) & The read is reduced to its 5' most base or 3' most base. Read summarization is then performed based on the single base position to which the read is reduced. By default, no read reduction will be performed.\\
\hline
diff --git a/src/HelperFunctions.c b/src/HelperFunctions.c
index f354645..39b9dd4 100644
--- a/src/HelperFunctions.c
+++ b/src/HelperFunctions.c
@@ -21,6 +21,27 @@
#include <string.h>
#include <assert.h>
+
+#ifdef MACOS
+
+#include <sys/types.h>
+#include <sys/socket.h>
+#include <sys/ioctl.h>
+#include <sys/sysctl.h>
+#include <net/if.h>
+#include <net/if_dl.h>
+#include <netinet/in.h>
+#include <arpa/inet.h>
+
+#else
+
+#include <sys/ioctl.h>
+#include <net/if.h>
+#include <unistd.h>
+#include <netinet/in.h>
+#endif
+
+
#include "subread.h"
#include "gene-algorithms.h"
#include "HelperFunctions.h"
@@ -740,3 +761,130 @@ int strcmp_number(char * s1, char * s2)
return ret;
}
}
+
+
+
+int mac_str(char * str_buff)
+{
+#ifdef FREEBSD
+ return 1;
+#else
+#ifdef MACOS
+ int mib[6], x1, ret = 1;
+ size_t len;
+ char *buf;
+ unsigned char *ptr;
+ struct if_msghdr *ifm;
+ struct sockaddr_dl *sdl;
+
+
+ for(x1 = 0 ; x1 < 40; x1++)
+ {
+ mib[0] = CTL_NET;
+ mib[1] = AF_ROUTE;
+ mib[2] = 0;
+ mib[3] = AF_LINK;
+ mib[4] = NET_RT_IFLIST;
+ mib[5] = x1;
+
+ if (sysctl(mib, 6, NULL, &len, NULL, 0) >=0) {
+ buf = malloc(len);
+ if (sysctl(mib, 6, buf, &len, NULL, 0) >=0) {
+
+ ifm = (struct if_msghdr *)buf;
+ sdl = (struct sockaddr_dl *)(ifm + 1);
+ ptr = (unsigned char *)LLADDR(sdl);
+
+ if(sdl -> sdl_nlen < 1) continue;
+
+ char * ifname = malloc(sdl -> sdl_nlen + 1);
+
+ memcpy(ifname, sdl -> sdl_data, sdl -> sdl_nlen);
+ ifname[sdl -> sdl_nlen] = 0;
+ if(ifname[0]!='e'){
+ free(ifname);
+ continue;
+ }
+ free(ifname);
+
+ sprintf(str_buff,"%02X%02X%02X%02X%02X%02X", *ptr, *(ptr+1), *(ptr+2),
+ *(ptr+3), *(ptr+4), *(ptr+5));
+ ret = 0;
+ break;
+ }
+ free(buf);
+ }
+ }
+ return ret;
+#else
+ struct ifreq ifr;
+ struct ifconf ifc;
+ char buf[1024];
+ int success = 0;
+
+ int sock = socket(AF_INET, SOCK_DGRAM, IPPROTO_IP);
+ if (sock == -1) { /* handle error*/ };
+
+ ifc.ifc_len = sizeof(buf);
+ ifc.ifc_buf = buf;
+ if (ioctl(sock, SIOCGIFCONF, &ifc) == -1) { /* handle error */ }
+
+ struct ifreq* it = ifc.ifc_req;
+ const struct ifreq* const end = it + (ifc.ifc_len / sizeof(struct ifreq));
+
+ for (; it != end; ++it) {
+ strcpy(ifr.ifr_name, it->ifr_name);
+ if (ioctl(sock, SIOCGIFFLAGS, &ifr) == 0) {
+ if (! (ifr.ifr_flags & IFF_LOOPBACK)) { // don't count loopback
+ if (ioctl(sock, SIOCGIFHWADDR, &ifr) == 0) {
+ success = 1;
+ break;
+ }
+ }
+ }
+ }
+
+ unsigned char mac_address[6];
+
+ if (success){
+ memcpy(mac_address, ifr.ifr_hwaddr.sa_data, 6);
+ int x1;
+ for(x1 = 0; x1 < 6; x1++){
+ sprintf(str_buff+2*x1, "%02X",mac_address[x1]);
+ }
+ return 0;
+ }
+ return 1;
+#endif
+#endif
+}
+
+int rand_str(char * str_buff){
+ int ret = 1;
+ FILE * fp = fopen("/dev/urandom","r");
+ if(fp){
+ int x1;
+ for(x1=0; x1<6; x1++){
+ sprintf(str_buff + 2*x1 , "%02X", fgetc(fp));
+ }
+ fclose(fp);
+ ret = 0;
+ }
+ return ret;
+}
+
+int mathrand_str(char * str_buff){
+ srand((int)(miltime()*100));
+ int x1;
+ for(x1 = 0; x1 < 6; x1++){
+ sprintf(str_buff+2*x1, "%02X", rand() & 0xff );
+ }
+ return 0;
+}
+
+int mac_or_rand_str(char * str_buff){
+ return mac_str(str_buff) && rand_str(str_buff) && mathrand_str(str_buff);
+}
+
+
+
diff --git a/src/HelperFunctions.h b/src/HelperFunctions.h
index defa736..2f40d21 100644
--- a/src/HelperFunctions.h
+++ b/src/HelperFunctions.h
@@ -68,4 +68,5 @@ int strcmp_number(char * s1, char * s2);
unsigned int reverse_cigar(unsigned int pos, char * cigar, char * new_cigar);
unsigned int find_left_end_cigar(unsigned int right_pos, char * cigar);
+int mac_or_rand_str(char * char_14);
#endif
diff --git a/src/Makefile b/src/Makefile
index 35f5a29..9215779 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -8,7 +8,7 @@ all:
@echo " For building subread in FreeBSD, please run ' gmake -f Makefile.FreeBSD '."
@echo " For building subread in Oracle Solaris or OpenSolaris, please run ' gmake -f Makefile.SunOS '."
@echo
- @echo " The default compiler is gcc; you may change the compailer by editing the makefile."
+ @echo " The default compiler is gcc; you may change it by editing the Makefiles for platforms."
@echo
@echo " The generated executables are saved to directory (PACKAGE_DIR)/bin"
@echo
diff --git a/src/Makefile.Linux b/src/Makefile.Linux
index abb246b..87ec810 100644
--- a/src/Makefile.Linux
+++ b/src/Makefile.Linux
@@ -1,10 +1,10 @@
-MACOS = -D MACOS
+#MACOS = -D MACOS
include makefile.version
-CCFLAGS = -mtune=core2 ${MACOS} -O9 -Wall -DMAKE_FOR_EXON -D MAKE_STANDALONE -D SUBREAD_VERSION=\"${SUBREAD_VERSION}\" -D_FILE_OFFSET_BITS=64
-#CCFLAGS = -D_FORTIFY_SOURCE=2 -mtune=core2 ${MACOS} -O2 -Wall -DMAKE_FOR_EXON -D MAKE_STANDALONE -D SUBREAD_VERSION=\"${SUBREAD_VERSION}\" # -DREPORT_ALL_THE_BEST
-LDFLAGS = ${STATIC_MAKE} -lpthread -lz -lm ${MACOS} -O9 -DMAKE_FOR_EXON -D MAKE_STANDALONE # -DREPORT_ALL_THE_BEST
+OPT_LEVEL = 9
+CCFLAGS = -mtune=core2 ${MACOS} -O${OPT_LEVEL} -Wall -DMAKE_FOR_EXON -D MAKE_STANDALONE -D SUBREAD_VERSION=\"${SUBREAD_VERSION}\" -D_FILE_OFFSET_BITS=64
+LDFLAGS = ${STATIC_MAKE} -lpthread -lz -lm ${MACOS} -O${OPT_LEVEL} -DMAKE_FOR_EXON -D MAKE_STANDALONE # -DREPORT_ALL_THE_BEST
CC = gcc ${CCFLAGS} -ggdb -fomit-frame-pointer -ffast-math -funroll-loops -mmmx -msse -msse2 -msse3 -fmessage-length=0
diff --git a/src/SNPCalling.c b/src/SNPCalling.c
index 7a77053..fbd8b03 100644
--- a/src/SNPCalling.c
+++ b/src/SNPCalling.c
@@ -32,6 +32,7 @@
#include <getopt.h>
#include "subread.h"
#include "hashtable.h"
+#include "HelperFunctions.h"
#include "gene-algorithms.h"
#include "SNPCalling.h"
#include "input-files.h"
@@ -506,8 +507,8 @@ void fishers_test_on_block(struct SNP_Calling_Parameters * parameters, float * s
flanking_matched = d;
}
- //SUBREADprintf("TEST: %u a,b,c,d=%d %d %d %d; FU=%d FM=%d; Goahead=%d; Tailleft=%d\n", i,a,b,c,d, flanking_unmatched, flanking_matched, go_ahead, left_tail);
float p_middle = fisher_exact_test(a, flanking_unmatched, c, flanking_matched);
+ //SUBREADprintf("TEST: %u a,b,c,d=%d %d %d %d; FU=%d FM=%d; Goahead=%d; Tailleft=%d; p=%G; p-cut=%G\n", i,a,b,c,d, flanking_unmatched, flanking_matched, go_ahead, left_tail, p_middle, p_cutoff);
if(all_result_needed || ( p_middle < p_cutoff && flanking_matched*20>(flanking_matched+ flanking_unmatched )*16))
snp_fisher_raw [i] = p_middle;
else snp_fisher_raw [i] = -999;
@@ -842,6 +843,8 @@ int process_snp_votes(FILE *out_fp, unsigned int offset , unsigned int reference
}
}
}
+
+ // SUBREADprintf("OUTTEXT: i=%d, all_reads=%d, Fisher-p=%f, SNPs=%d, POI-MM=%d\n", i, all_reads, snp_fisher_raw[i], snps, POI_MM);
if(snps && POI_MM *1. / all_reads >= parameters->supporting_read_rate )
{
if(snp_fisher_raw[i] >= 0. || parameters -> fisher_exact_testlen < 1)
@@ -1008,8 +1011,8 @@ int run_chromosome_search(FILE *in_fp, FILE * out_fp, char * chro_name , char *
if( (*task_no) % all_threads == thread_no && all_offset <= chro_len)
{
- //#warning "=== ONLY TEST ONE BLOCK ==="
- //if(strcmp(chro_name,"chr7")==0 && all_offset == 45000000){
+ //#warning "=== ONLY TEST ONE BLOCK , USE 'if(1)' IN RELEASE ==="
+ //if(strcmp(chro_name,"chr7")==0 && all_offset == 60000000){
if(1){
process_snp_votes(out_fp, all_offset, offset, referenced_base, chro_name , temp_prefix, parameters);
print_in_box(89,0,0,"processed block %c[36m%s@%d%c[0m by thread %d/%d [block number=%d/%d]", CHAR_ESC, chro_name, all_offset, CHAR_ESC , thread_no+1, all_threads, 1+(*task_no)-parameters->empty_blocks, parameters->all_blocks);
@@ -1018,7 +1021,7 @@ int run_chromosome_search(FILE *in_fp, FILE * out_fp, char * chro_name , char *
else if((*task_no) % all_threads == thread_no)
{
print_in_box(80,0,0,"Ignored in: %s@%d by thr %d/%d [tid=%d]\n", chro_name, all_offset, thread_no, all_threads, *task_no);
- SUBREADprintf("LEN %u > %u\n", all_offset, chro_len);
+ // SUBREADprintf("LEN %u > %u\n", all_offset, chro_len);
}
offset = 0;
(*task_no)++;
@@ -1382,11 +1385,14 @@ int SNP_calling(char * in_SAM_file, char * out_BED_file, char * in_FASTA_file, c
HashTableSetKeyComparisonFunction(parameters-> cigar_event_table, my_strcmp);
memcpy(rand48_seed, &start_time, 6);
- srand(time(NULL));
if(temp_location)
strcpy(temp_file_prefix, temp_location);
- else
- sprintf(temp_file_prefix, "./temp-snps-%06u-%08X-", getpid(), rand());
+ else{
+ char mac_rand[13];
+ mac_or_rand_str(mac_rand);
+
+ sprintf(temp_file_prefix, "./temp-snps-%06u-%s-", getpid(), mac_rand);
+ }
_EXSNP_SNP_delete_temp_prefix = temp_file_prefix;
print_in_box(89,0,0,"Split %s file into %c[36m%s*%c[0m ..." , parameters -> is_BAM_file_input?"BAM":"SAM" , CHAR_ESC, temp_file_prefix, CHAR_ESC);
diff --git a/src/core-bigtable.c b/src/core-bigtable.c
index 4cfc1a2..52cd5a8 100644
--- a/src/core-bigtable.c
+++ b/src/core-bigtable.c
@@ -86,21 +86,23 @@ int init_bigtable_results(global_context_t * global_context, int is_rewinding)
global_context -> bigtable_chunked_fragments = global_context -> config.reads_per_chunk+1;
global_context -> bigtable_cache_size = global_context -> bigtable_chunked_fragments * (1+global_context -> input_reads.is_paired_end_reads);
} else {
- global_context -> bigtable_chunked_fragments = 300000 - 260000;
- global_context -> bigtable_cache_size = global_context -> config.all_threads * global_context -> bigtable_chunked_fragments * (1+global_context -> input_reads.is_paired_end_reads);
+ global_context -> bigtable_chunked_fragments = 110llu*1024*1024;
+ global_context -> bigtable_cache_size = (1+0*global_context -> config.all_threads) * global_context -> bigtable_chunked_fragments * (1+global_context -> input_reads.is_paired_end_reads);
}
//SUBREADprintf("reads_per_chunk = %u ; cached_single_reads = %u ; size of each read = %d + %d\n", global_context -> config.reads_per_chunk, global_context -> bigtable_cache_size, sizeof(mapping_result_t) , sizeof(subjunc_result_t));
- if(!is_rewinding)
+ if(!is_rewinding){
global_context -> bigtable_cache = malloc(sizeof(bigtable_cached_result_t) * global_context -> bigtable_cache_size);
-
+ }
int xk1;
for(xk1 = 0; xk1 < global_context -> bigtable_cache_size; xk1++){
- if(!is_rewinding)
+ if(!is_rewinding){
global_context -> bigtable_cache [xk1].alignment_res = malloc(sizeof(mapping_result_t) * global_context -> config.multi_best_reads);
+ assert(global_context -> bigtable_cache [xk1].alignment_res!=NULL);
+ }
if(global_context -> config.use_memory_buffer)
{
@@ -134,6 +136,7 @@ int init_bigtable_results(global_context_t * global_context, int is_rewinding)
return 0;
}
+
mapping_result_t * _global_retrieve_alignment_ptr(global_context_t * global_context, subread_read_number_t pair_number, int is_second_read, int best_read_id){
mapping_result_t * ret;
bigtable_retrieve_result(global_context, NULL, pair_number, best_read_id, is_second_read, &ret, NULL);
@@ -191,8 +194,13 @@ bigtable_cached_result_t * bigtable_retrieve_cache(global_context_t * global_con
long long load_start_pair_no = inner_pair_number - inner_pair_number % global_context -> bigtable_chunked_fragments;
if(global_context -> bigtable_cache_file_fp){
+
+
+ //SUBREADprintf("MARK_OCCPY=%lld BY THREAD %d\n", pair_number, thread_context ? thread_context -> thread_id : -1);
+
if(global_context -> bigtable_cache_file_loaded_fragments_begin == -1 || inner_pair_number >= global_context -> bigtable_cache_file_loaded_fragments_begin + global_context -> bigtable_chunked_fragments || inner_pair_number < global_context -> bigtable_cache_file_loaded_fragments_begin)
{
+ SUBREADprintf("THREAD # %d WAITING FOR %llu for RETRIEVE %llu\n", thread_context? thread_context -> thread_id:-1, global_context -> bigtable_cache_file_loaded_fragments_begin, pair_number);
wait_occupied(global_context, global_context -> bigtable_cache_file_loaded_fragments_begin);
}
@@ -238,7 +246,7 @@ bigtable_cached_result_t * bigtable_retrieve_cache(global_context_t * global_con
if(0 && xk1 < 10)
{
//SUBREADprintf("CACHEP_211: %p (%d from %llu)\n", current_cache, xk1, pair_number);
- SUBREADprintf("NENSET_211: %p\n", current_cache -> alignment_res);
+ SUBREADprintf("NENSET_211: %p + %d\n", current_cache -> alignment_res, xk1 * (1+global_context -> input_reads.is_paired_end_reads) + xk2);
}
memset( current_cache -> alignment_res , 0, sizeof(mapping_result_t) * global_context -> config.multi_best_reads);
@@ -259,8 +267,9 @@ bigtable_cached_result_t * bigtable_retrieve_cache(global_context_t * global_con
bigtable_unlock(global_context);
}
- if(global_context -> bigtable_cache_file_fp)
+ if(global_context -> bigtable_cache_file_fp){
global_context -> bigtable_cache[inner_pair_number - load_start_pair_no].status = CACHE_STATUS_OCCUPIED;
+ }
bigtable_cached_result_t * ret_cache = global_context -> bigtable_cache + (inner_pair_number - load_start_pair_no)* (1+global_context -> input_reads.is_paired_end_reads) + is_second_read;
return ret_cache;
@@ -285,15 +294,17 @@ unsigned short * _global_retrieve_big_margin_ptr(global_context_t * global_conte
}
-// do nothing : the data is automatically saved into the temporary file when using mmap.
void bigtable_release_result(global_context_t * global_context , thread_context_t * thread_context , subread_read_number_t pair_number, int commit_change){
- long long inner_pair_number = get_inner_pair(global_context, pair_number);
- long long load_start_pair_no = inner_pair_number - inner_pair_number % global_context -> bigtable_chunked_fragments;
- if(global_context -> bigtable_cache_file_fp)
+ if(global_context -> bigtable_cache_file_fp){
+ long long inner_pair_number = get_inner_pair(global_context, pair_number);
+ long long load_start_pair_no = inner_pair_number - inner_pair_number % global_context -> bigtable_chunked_fragments;
+ //SUBREADprintf("REAL RELEASE:%lld\n", inner_pair_number - load_start_pair_no);
global_context -> bigtable_cache[inner_pair_number - load_start_pair_no].status = CACHE_STATUS_RELEASED;
+ //SUBREADprintf("OCCUPYXX0 = %d\n", global_context -> bigtable_cache[0].status);
- if(commit_change){
- global_context -> bigtable_dirty_data=1;
+ if(commit_change){
+ global_context -> bigtable_dirty_data=1;
+ }
}
}
@@ -479,8 +490,11 @@ void fraglist_destroy(fragment_list_t * list){
void fraglist_append(fragment_list_t * list, subread_read_number_t fragment_number){
if(list -> fragments >= list -> capacity){
+ //#warning "============== COMMENT DEBUG INFO ====================="
+ //SUBREADprintf("REALLOC_PRT_IN = %d , %p\n",list -> capacity , list -> fragment_numbers );
list -> capacity = max(list -> capacity + 5, list -> capacity * 1.3);
list -> fragment_numbers = realloc(list -> fragment_numbers, sizeof(subread_read_number_t) * list -> capacity);
+ // SUBREADprintf("REALLOC_PRT_OUT = %d , %p\n",list -> capacity , list -> fragment_numbers );
}
list -> fragment_numbers[ list -> fragments ++ ] = fragment_number;
diff --git a/src/core-indel.c b/src/core-indel.c
index 4fae9a8..9a0fe5a 100644
--- a/src/core-indel.c
+++ b/src/core-indel.c
@@ -25,6 +25,7 @@
#include <assert.h>
#include <unistd.h>
#include "hashtable.h"
+#include "HelperFunctions.h"
#include "gene-value-index.h"
#include "gene-algorithms.h"
#include "input-files.h"
@@ -1981,7 +1982,7 @@ int write_local_reassembly(global_context_t *global_context, HashTable *pileup_f
- if(0 == locate_gene_position_max(anchor_pos, &global_context -> chromosome_table, &chro_name, &chro_offset, read_len))
+ if(0 == locate_gene_position_max(anchor_pos, &global_context -> chromosome_table, &chro_name, &chro_offset, NULL, NULL, read_len))
{
char temp_file_name[MAX_FILE_NAME_LENGTH];
int close_now = 0;
@@ -3654,8 +3655,8 @@ int finalise_pileup_file_by_voting(global_context_t * global_context , char * te
char * chro_begin;
int chro_offset_start = 0;
int chro_offset_end = 0;
- locate_gene_position_max(contig_start_pos + head_removed_bases ,& global_context -> chromosome_table, &chro_begin, &chro_offset_start, 0);
- locate_gene_position_max(contig_end_pos - tail_removed_bases ,& global_context -> chromosome_table, &chro_begin, &chro_offset_end, 0);
+ locate_gene_position_max(contig_start_pos + head_removed_bases ,& global_context -> chromosome_table, &chro_begin, &chro_offset_start, NULL, NULL, 0);
+ locate_gene_position_max(contig_end_pos - tail_removed_bases ,& global_context -> chromosome_table, &chro_begin, &chro_offset_end, NULL, NULL, 0);
if(full_rebuilt_window_size - read_position_cursor - tail_removed_bases)
fprintf(global_context -> long_insertion_FASTA_fp, ">%s-%u-%u-%s%dM\n", chro_name, chro_offset_start, chro_offset_end - 1, contig_CIGAR, full_rebuilt_window_size - read_position_cursor - tail_removed_bases);
else
@@ -4003,9 +4004,6 @@ void init_global_context(global_context_t * context)
context->config.show_soft_cliping = 1 ;
context->config.big_margin_record_size = 9;
- //#warning "====== FOR HIGH JUNCTION ACCURACT ==========="
- context->config.big_margin_record_size = 15;
-
context->config.read_group_id[0] = 0;
context->config.read_group_txt[0] = 0;
context->config.first_read_file[0] = 0;
@@ -4033,7 +4031,10 @@ void init_global_context(global_context_t * context)
memcpy(seed_rand, &double_time, 2*sizeof(int));
srand(seed_rand[0]^seed_rand[1]);
- sprintf(context->config.temp_file_prefix, "./core-temp-sum-%06u-%06u", getpid(),rand());
+ char mac_rand[13];
+ mac_or_rand_str(mac_rand);
+
+ sprintf(context->config.temp_file_prefix, "./core-temp-sum-%06u-%s", getpid(), mac_rand );
_COREMAIN_delete_temp_prefix = context->config.temp_file_prefix;
@@ -4060,6 +4061,7 @@ void init_global_context(global_context_t * context)
int CORE_DPALIGN_MATCH_SCORE = 2;
int CORE_DPALIGN_MISMATCH_PENALTY = 0;
+
int core_dynamic_align(global_context_t * global_context, thread_context_t * thread_context, char * read, int read_len, unsigned int begin_position, char * movement_buffer, int expected_offset, char * read_name)
// read must be converted to the positive strand.
// movement buffer: 0:match, 1: read-insert, 2: gene-insert, 3:mismatch
@@ -4081,6 +4083,9 @@ int core_dynamic_align(global_context_t * global_context, thread_context_t * thr
gene_value_index_t * current_value_index = thread_context?thread_context->current_value_index:global_context->current_value_index;
indel_context_t * indel_context = (indel_context_t*)global_context -> module_contexts[MODULE_INDEL_ID];
+
+ //unsigned long long table_ptr = (unsigned long long) indel_context -> dynamic_align_table;
+
short ** table = indel_context -> dynamic_align_table;
char ** table_mask = indel_context -> dynamic_align_table_mask;
if(thread_context)
diff --git a/src/core-interface-aligner.c b/src/core-interface-aligner.c
index 1093f93..e870452 100644
--- a/src/core-interface-aligner.c
+++ b/src/core-interface-aligner.c
@@ -467,6 +467,8 @@ int parse_opts_aligner(int argc , char ** argv, global_context_t * global_contex
else if(strcmp("SVdetection", long_options[option_index].name)==0)
{
global_context -> config.do_structural_variance_detection = 1;
+ global_context -> config.use_memory_buffer = 1;
+ global_context -> config.reads_per_chunk = 600llu*1024*1024;
}
else if(strcmp("maxRealignLocations", long_options[option_index].name)==0)
{
diff --git a/src/core-interface-subjunc.c b/src/core-interface-subjunc.c
index e0b316c..a9d43ef 100644
--- a/src/core-interface-subjunc.c
+++ b/src/core-interface-subjunc.c
@@ -60,7 +60,6 @@ static struct option long_options[] =
{"maxRealignLocations", required_argument, 0, 0},
{"minVoteCutoff", required_argument, 0, 0},
{"minMappedFraction", required_argument, 0, 0},
- {"disableBigMargin", no_argument, 0, 0},
{"complexIndels", no_argument, 0, 0},
{0, 0, 0, 0}
};
@@ -482,6 +481,8 @@ int parse_opts_subjunc(int argc , char ** argv, global_context_t * global_contex
else if(strcmp("SVdetection", long_options[option_index].name)==0)
{
global_context -> config.do_structural_variance_detection = 1;
+ global_context -> config.use_memory_buffer = 1;
+ global_context -> config.reads_per_chunk = 300llu*1024*1024;
}
else if(strcmp("minDistanceBetweenVariants", long_options[option_index].name)==0)
{
@@ -540,6 +541,23 @@ int parse_opts_subjunc(int argc , char ** argv, global_context_t * global_contex
if(global_context->config.is_SAM_file_input) global_context->config.phred_score_format = FASTQ_PHRED33;
global_context->config.more_accurate_fusions = global_context->config.more_accurate_fusions && global_context->config.do_fusion_detection;
+ if(global_context->config.more_accurate_fusions)
+ {
+ global_context->config.high_quality_base_threshold = 999999;
+ //#warning "============ REMOVE THE NEXT LINE ======================"
+ global_context->config.show_soft_cliping = 1;
+ //#warning "============ REMOVE ' + 3' FROM NEXT LINE =============="
+ global_context->config.max_mismatch_junction_reads = 0 + 3;
+
+ //#warning "============ REMOVE ' - 1' FROM NEXT LINE =============="
+ global_context->config.do_big_margin_filtering_for_junctions = 1 - 1;
+ global_context->config.total_subreads = 28;
+
+ //#warning "============ REMOVE THE NEXT LINE BEFORE RELEASE ==============="
+ //global_context->config.multi_best_reads = 1;
+ }
+
+
return 0;
}
diff --git a/src/core-junction.c b/src/core-junction.c
index 7ff76ba..3e94dbc 100644
--- a/src/core-junction.c
+++ b/src/core-junction.c
@@ -34,6 +34,11 @@
#define TTTSNAME "V0112_0155:7:1308:1308:136442"
+#define CLUSTER_ALIGNMENT_DONOR_R1_MAPPED 2
+#define CLUSTER_ALIGNMENT_DONOR_R2_MAPPED 4
+#define CLUSTER_ALIGNMENT_DONOR_NEGATIVE_STRAND 1
+
+
unsigned int abs32uint(unsigned int x){
if(x > 0x7fffffff) x = (0xffffffff - x) + 1;
return x;
@@ -103,6 +108,9 @@ typedef struct{
void search_events_to_front(global_context_t * global_context, thread_context_t * thread_context, explain_context_t * explain_context, char * read_text , char * qual_text, unsigned int read_head_abs_offset, short remainder_len, short sofar_matched, int suggested_movement, int do_not_jump)
{
short tested_read_pos;
+ // #warning "SUBREAD_151 REMOVE THIS ASSERTION! "
+ // if(remainder_len >= 102)SUBREADprintf("FATAL:%d\n", remainder_len );
+ // assert(remainder_len < 102);
HashTable * event_table = NULL;
chromosome_event_t * event_space = NULL;
@@ -223,6 +231,8 @@ void search_events_to_front(global_context_t * global_context, thread_context_t
short new_remainder_len = remainder_len - tested_read_pos + min(0, tested_event->indel_length) - tested_event -> indel_at_junction;
+ // #warning "SUBREAD_151 REMOVE THIS ASSERTION! "
+ // assert(new_remainder_len < 102);
if(new_remainder_len>0)
{
@@ -263,6 +273,8 @@ void search_events_to_front(global_context_t * global_context, thread_context_t
//printf("JUMP_IN: %u ; STRAND=%c ; REMENDER=%d ; 0=%d 0=%d\n", new_read_head_abs_offset, tested_event -> is_strand_jumped?'X':'=', new_remainder_len, tested_event -> indel_length, tested_event -> indel_at_junction);
//}
+ // #warning "SUBREAD_151 REMOVE THIS ASSERTION! "
+ // assert(new_remainder_len < 102);
//printf("SUGGEST_NEXT = %d (! %d)\n", tested_event -> connected_next_event_distance, tested_event -> connected_previous_event_distance);
search_events_to_front(global_context, thread_context, explain_context, read_text + tested_event -> indel_at_junction + tested_read_pos - min(0, tested_event->indel_length), qual_text + tested_read_pos - min(0, tested_event->indel_length), new_read_head_abs_offset, new_remainder_len, sofar_matched + matched_bases_to_site - jump_penalty, tested_event -> connected_next_event_distance, 0);
explain_context -> tmp_search_sections --;
@@ -279,6 +291,8 @@ void search_events_to_front(global_context_t * global_context, thread_context_t
is_junction_scanned = max(is_junction_scanned, this_round_junction_scanned);
}
+ //#warning "SUBREAD_151 REMOVE THE ASSERT! "
+ //assert( remainder_len< 102 );
int whole_section_matched = match_chro(read_text , value_index, explain_context -> current_is_strand_jumped?read_head_abs_offset - remainder_len +1:read_head_abs_offset, remainder_len , explain_context -> current_is_strand_jumped, global_context -> config.space_type);
explain_context -> tmp_total_matched_bases = whole_section_matched + sofar_matched ;
@@ -387,21 +401,27 @@ void new_explain_try_replace(global_context_t* global_context, thread_context_t
if(search_to_back){
explain_context -> all_back_alignments = 1;
explain_context -> result_back_junction_numbers[0] = explain_context -> tmp_search_sections +1;
+ // checked: memory boundary
memcpy(explain_context -> result_back_junctions[0], explain_context -> tmp_search_junctions , sizeof(perfect_section_in_read_t) * (explain_context -> tmp_search_sections +1));
}else{
explain_context -> all_front_alignments = 1;
explain_context -> result_front_junction_numbers[0] = explain_context -> tmp_search_sections +1;
+ // checked: memory boundary
memcpy(explain_context -> result_front_junctions[0], explain_context -> tmp_search_junctions , sizeof(perfect_section_in_read_t) * (explain_context -> tmp_search_sections +1));
}
}else if(is_same_best){
if(search_to_back && explain_context -> all_back_alignments < MAX_ALIGNMENT_PER_ANCHOR){
explain_context -> result_back_junction_numbers[explain_context -> all_back_alignments] = explain_context -> tmp_search_sections +1;
+
+ // checked: memory boundary
memcpy(explain_context -> result_back_junctions[explain_context -> all_back_alignments], explain_context -> tmp_search_junctions , sizeof(perfect_section_in_read_t) * (explain_context -> tmp_search_sections +1));
explain_context -> all_back_alignments ++;
}else if((!search_to_back) && explain_context -> all_front_alignments < MAX_ALIGNMENT_PER_ANCHOR){
explain_context -> result_front_junction_numbers[explain_context -> all_front_alignments] = explain_context -> tmp_search_sections +1;
+
+ // checked: memory boundary
memcpy(explain_context -> result_front_junctions[explain_context -> all_front_alignments], explain_context -> tmp_search_junctions , sizeof(perfect_section_in_read_t) * (explain_context -> tmp_search_sections +1));
explain_context -> all_front_alignments ++;
}
@@ -670,103 +690,6 @@ void insert_big_margin_record(global_context_t * global_context , unsigned short
}
}
-// This function try to add a new anchor into the list or replace an existing anchor by moving done the following anchors.
-// It is only invoked in the first step: select the best anchors. No minor half is considered at all.
-// It also makes if the current result is a tie score: if the last and current Vote+Coverage+Hamming+Qual are equal
-void do_append_inner(global_context_t * global_context, thread_context_t * thread_context, subread_read_number_t pair_number, int * used_anchors, int total_anchors, select_junction_record_t * anchor_list, gene_vote_number_t Vote_major, int coverage_major_start, int coverage_major_end, int hamming_major, int quality_major, unsigned int pos_major, int flags, int read_len, gene_vote_number_t * indel_recorder)
-{
- int xx;
- int replace_index = -1;
- int i_am_break_even = 0;
- if(0<*used_anchors)
- {
- for(xx=0; xx< *used_anchors;xx++){
- select_junction_record_t * tanchor = anchor_list + xx;
-
- if(Vote_major >tanchor -> piece_main_votes ||
- (Vote_major ==tanchor -> piece_main_votes && coverage_major_end-coverage_major_start > tanchor -> piece_main_coverage_end-tanchor -> piece_main_coverage_start) ||
- (Vote_major ==tanchor -> piece_main_votes && coverage_major_end-coverage_major_start ==tanchor -> piece_main_coverage_end-tanchor -> piece_main_coverage_start && hamming_major > tanchor -> piece_main_hamming_match) ||
- (Vote_major ==tanchor -> piece_main_votes && coverage_major_end-coverage_major_start ==tanchor -> piece_main_coverage_end-tanchor -> piece_main_coverage_start && hamming_major ==tanchor -> piece_main_hamming_match && quality_major >= tanchor -> piece_main_read_quality))
- {
- if((Vote_major ==tanchor -> piece_main_votes && coverage_major_end-coverage_major_start ==tanchor -> piece_main_coverage_end-tanchor -> piece_main_coverage_start && hamming_major ==tanchor -> piece_main_hamming_match && quality_major == tanchor -> piece_main_read_quality))
- {
- // a tie
- if(xx < total_anchors - 1)
- replace_index = xx;
-
- if(xx == 0)// the BEST anchor is a tie
- {
- tanchor -> is_break_even = 1;
-
- int yy;
- for(yy = 1; yy < *used_anchors; yy++)
- {
- select_junction_record_t * canchor = anchor_list + yy;
- if((Vote_major ==canchor -> piece_main_votes && coverage_major_end-coverage_major_start ==canchor -> piece_main_coverage_end-canchor -> piece_main_coverage_start && hamming_major ==canchor -> piece_main_hamming_match && quality_major == canchor -> piece_main_read_quality))
- canchor -> is_break_even = 1;
- }
- i_am_break_even = 1;
- }
-
- break;
- }
- else
- {
- // the current XX-th item is move down.
- replace_index = xx;
- if(xx == 0) // the BEST anchor is clearly not a tie
- {
- int yy;
- for(yy = 0; yy < *used_anchors; yy++)
- {
- select_junction_record_t * canchor = anchor_list + yy;
- canchor -> is_break_even = 0;
- }
- }
- break;
- }
- }
-
- if(replace_index < 0 && (*used_anchors) < total_anchors ) replace_index = (*used_anchors);
- }
- }else replace_index = 0;
-
-
- if(replace_index >= 0){
- for(xx = (* used_anchors) - 1; xx >= replace_index ; xx--)
- {
- if(xx < total_anchors - 1)
- memcpy(anchor_list + xx+1, anchor_list+xx, sizeof( select_junction_record_t ));
- }
-
- int major_indels = 0;
-
- if(read_len > EXON_LONG_READ_LENGTH){
- int kx1;
- for(kx1=0; kx1<MAX_INDEL_SECTIONS; kx1++)
- {
- if(!indel_recorder[kx1*3]) break;
- major_indels += indel_recorder[kx1*3+2];
- }
- }
-
- select_junction_record_t * nanchor = anchor_list + replace_index;
- memset(nanchor , 0 , sizeof( select_junction_record_t ));
- nanchor -> is_break_even = i_am_break_even;
- nanchor -> piece_main_votes = Vote_major;
- nanchor -> piece_main_coverage_start = coverage_major_start;
- nanchor -> piece_main_coverage_end = coverage_major_end;
- nanchor -> piece_main_hamming_match = hamming_major;
- nanchor -> piece_main_read_quality = quality_major;
- nanchor -> piece_main_abs_offset = pos_major;
- nanchor -> piece_main_masks = flags;
- nanchor -> piece_main_indels = major_indels;
- nanchor -> piece_main_indel_record = indel_recorder;
-
- if( * used_anchors < total_anchors) (*used_anchors) ++;
- }
-}
-
int is_PE_distance(global_context_t * global_context, unsigned int pos1, unsigned int pos2, int rlen1, int rlen2, int is_negative_R1, int is_negative_R2)
{
long long int dist = pos2;
@@ -1038,7 +961,7 @@ void copy_vote_to_alignment_res(global_context_t * global_context, thread_contex
if(global_context->config.do_fusion_detection){
// function test_small_minor_votes returns 1 if the vote number is not significantly
// higher than the vote numbers in the vote list.
- //#warning "=========== THE TWO LINES SHOULD BE UNCOMMENTED IN RELEASED VERSION ==== WE COMMENT IT FOR A BETTER FUSION SENSITIVITY BUT ONLY FOR TEST ==================="
+ //#warning "SUBREAD_151 =========== THE TWO LINES SHOULD BE UNCOMMENTED IN RELEASED VERSION ==== WE COMMENT IT FOR A BETTER FUSION SENSITIVITY BUT ONLY FOR TEST ==================="
if(1){
int small_minor_bigmargin = test_small_minor_votes(global_context , i, j, vote_i, vote_j, current_vote, curr_read_len);
if(small_minor_bigmargin) continue;
@@ -1206,14 +1129,15 @@ void copy_vote_to_alignment_res(global_context_t * global_context, thread_contex
if(replace_minor>0) replace_minor += current_vote -> votes[i][j] * 100000;
- //SUBREADprintf("NOJUMP_DONORs=%d LOC=%u\n", replace_minor , current_vote -> pos[i][j]);
+ if(0 && ( FIXLENstrcmp("R006232475", read_name) == 0 ) )
+ SUBREADprintf("NOJUMP_DONORs=%d LOC=%u\n", replace_minor , current_vote -> pos[i][j]);
if(global_context -> config.do_fusion_detection && !(current_vote -> masks[vote_i][vote_j] & IS_NEGATIVE_STRAND))
// changed back.
reverse_read(curr_read_text, curr_read_len, global_context -> config.space_type);
}
}
- if(0 && ( FIXLENstrcmp("R001968841", read_name) == 0 || FIXLENstrcmp("XXR002430582", read_name) == 0 ) )
+ if(0 && ( FIXLENstrcmp("R006232475", read_name) == 0 ) )
{
char posout[100];
absoffset_to_posstr(global_context, current_vote -> pos[i][j], posout);
@@ -1282,6 +1206,870 @@ void simple_PE_and_same_chro(global_context_t * global_context , simple_mapping_
test_PE_and_same_chro(global_context, r1 -> mapping_position, r2 -> mapping_position, is_PE_distance, is_same_chromosome, rlen1, rlen2);
}
+
+#define MAX_CLUSTER_ELEMENTS 7
+
+struct cluster_element{
+ unsigned int initial_position;
+ char cluster_members;
+ char from_second_read[MAX_CLUSTER_ELEMENTS];
+ int i_list[MAX_CLUSTER_ELEMENTS], j_list[MAX_CLUSTER_ELEMENTS];
+};
+
+int add_cluster_member(struct cluster_element * cl , int i, int j, int is_second_read){
+ if(cl->cluster_members < MAX_CLUSTER_ELEMENTS){
+ cl->i_list[(int)cl->cluster_members] = i;
+ cl->j_list[(int)cl->cluster_members] = j;
+ cl->from_second_read[(int)cl->cluster_members] = is_second_read;
+ cl->cluster_members++;
+ }
+ return cl->cluster_members;
+}
+
+int is_same_cluster( global_context_t * global_context, struct cluster_element * cl , unsigned int pos){
+ long long int test_pos = pos;
+ test_pos -= cl -> initial_position;
+ if(abs(test_pos) < global_context -> config.maximum_intron_length)
+ return 1;
+ return 0;
+}
+
+int process_voting_junction_PE_topK(global_context_t * global_context, thread_context_t * thread_context, subread_read_number_t pair_number, gene_vote_t * vote_1, gene_vote_t * vote_2, char * read_name_1, char * read_name_2, char * read_text_1, char * read_text_2, int read_len_1, int read_len_2, int is_negative_strand, gene_vote_number_t v1_all_subreads, gene_vote_number_t v2_all_subreads);
+int align_cluster(global_context_t * global_context, thread_context_t * thread_context, struct cluster_element * this_cluster, char * read_name_1, char * read_name_2, char * read_text_1, char * read_text_2, int read_len_1, int read_len_2, int is_negative_strand, gene_vote_t * vote_1, gene_vote_t * vote_2, int * this_score, int * ii_path, int * jj_path, int * masks, int * path_len, int * R1R2_mapped);
+
+void simple_copy_vote_to_result( mapping_result_t * align_res, gene_vote_t * current_vote, int vote_i, int vote_j, int used_subreads_in_vote, int noninformative_subreads_in_vote, int score){
+ align_res -> selected_position = current_vote -> pos[vote_i][vote_j];
+ align_res -> selected_votes = score;
+ align_res -> indels_in_confident_coverage = indel_recorder_copy(align_res -> selected_indel_record, current_vote -> indel_recorder[vote_i][vote_j]);
+ align_res -> confident_coverage_end = current_vote -> coverage_end[vote_i][vote_j];
+ align_res -> confident_coverage_start = current_vote -> coverage_start[vote_i][vote_j];
+ align_res -> result_flags = (current_vote -> masks[vote_i][vote_j] & IS_NEGATIVE_STRAND)?(CORE_IS_NEGATIVE_STRAND):0;
+ align_res -> used_subreads_in_vote = used_subreads_in_vote;
+ align_res -> noninformative_subreads_in_vote = noninformative_subreads_in_vote;
+}
+
+int process_voting_junction_PE_juncs( global_context_t * global_context, thread_context_t * thread_context, subread_read_number_t pair_number, gene_vote_t * vote_1, gene_vote_t * vote_2, char * read_name_1, char * read_name_2, char * read_text_1, char * read_text_2, int read_len_1, int read_len_2, int is_negative_strand, gene_vote_number_t v1_all_subreads, gene_vote_number_t v2_all_subreads ){
+ int current_cluster_number = 0,max_clusters = global_context -> config.max_vote_simples * 2;
+ int i,j, is_second_read, tested_votes, x1;
+
+ struct cluster_element * cluster_buffer = malloc(max_clusters * sizeof(struct cluster_element));
+ int max_cluster_size_r1 = 0, max_cluster_size_r2 = 0;
+
+ for( tested_votes = max(vote_1 -> max_vote, vote_2 -> max_vote); tested_votes > 0; tested_votes--) {
+ for(is_second_read = 0 ; is_second_read < 1 + global_context -> input_reads.is_paired_end_reads; is_second_read ++) {
+ gene_vote_t * current_vote = is_second_read?vote_2:vote_1;
+ int * max_cluster_size = is_second_read?(&max_cluster_size_r2):(&max_cluster_size_r1);
+ for (i=0; i<GENE_VOTE_TABLE_SIZE; i++) {
+ for (j=0; j< current_vote->items[i]; j++) {
+ if(current_vote->votes[i][j]!=tested_votes) continue;
+ int is_added = 0;
+
+ for(x1 = 0; x1 < current_cluster_number ; x1++){
+ if(is_same_cluster(global_context, cluster_buffer+x1, current_vote->pos[i][j])){
+ int new_size =add_cluster_member(cluster_buffer+x1, i, j, is_second_read);
+ (*max_cluster_size) = max(*max_cluster_size, new_size);
+ is_added = 1;
+ }
+ if(is_added)break;
+ }
+ if(current_cluster_number < max_clusters && !is_added){
+ cluster_buffer[current_cluster_number].initial_position = current_vote->pos[i][j];
+ cluster_buffer[current_cluster_number].i_list[0] = i;
+ cluster_buffer[current_cluster_number].j_list[0] = j;
+ cluster_buffer[current_cluster_number].from_second_read[0] = is_second_read;
+ cluster_buffer[current_cluster_number].cluster_members = 1;
+ current_cluster_number++;
+ }
+ }
+ }
+ }
+ }
+
+ if(1 || max_cluster_size_r1 == 3 || max_cluster_size_r2 == 3 ) // if there are 3-section clusters then parse it, else go to the regular procedure. There is a upper-limit to the sections to avoid fragile mapping.
+ {
+ for(x1 = 0 ; x1 < current_cluster_number ; x1++){
+ int this_score = -1, path_len = -1, R1R2_mapped = 0;
+ int this_ii_path[ MAX_CLUSTER_ELEMENTS ], this_jj_path[ MAX_CLUSTER_ELEMENTS ], this_masks [ MAX_CLUSTER_ELEMENTS ];
+ align_cluster(global_context, thread_context, cluster_buffer + x1, read_name_1, read_name_2, read_text_1, read_text_2, read_len_1, read_len_2, is_negative_strand, vote_1, vote_2, &this_score, this_ii_path, this_jj_path, this_masks, &path_len, &R1R2_mapped);
+
+ if(0 && FIXLENstrcmp("V0112_0155:7:1101:7309:2770", read_name_1)==0)
+ SUBREADprintf("REAE_TEST : R12MAP=%d, PATHLEN=%d, SCORE=%d\n", R1R2_mapped, path_len, this_score);
+
+ if(this_score > 0){
+ if(( R1R2_mapped & CLUSTER_ALIGNMENT_DONOR_R1_MAPPED) && ( R1R2_mapped & CLUSTER_ALIGNMENT_DONOR_R2_MAPPED)){
+ for(i = 0; i < global_context -> config.multi_best_reads; i++) {
+ mapping_result_t * old_result_R1 = _global_retrieve_alignment_ptr(global_context, pair_number, 0, i);
+ mapping_result_t * old_result_R2 = _global_retrieve_alignment_ptr(global_context, pair_number, 1, i);
+ short old_score_R1 = old_result_R1 -> selected_votes;
+ short old_score_R2 = old_result_R2 -> selected_votes;
+
+ if( old_score_R1 < this_score || old_score_R2 < this_score ){
+
+ for(j = global_context -> config.multi_best_reads - 2; j>=i; j--){
+ mapping_result_t * shifted_result_R1 = _global_retrieve_alignment_ptr(global_context, pair_number, 0, j);
+ mapping_result_t * shifted_result_R2 = _global_retrieve_alignment_ptr(global_context, pair_number, 1, j);
+ mapping_result_t * target_result_R1 = _global_retrieve_alignment_ptr(global_context, pair_number, 0, j+1);
+ mapping_result_t * target_result_R2 = _global_retrieve_alignment_ptr(global_context, pair_number, 1, j+1);
+ memcpy( target_result_R1, shifted_result_R1 , sizeof(mapping_result_t));
+ memcpy( target_result_R2, shifted_result_R2 , sizeof(mapping_result_t) );
+
+ }
+
+ int best_R1_i = -1, best_R1_j = - 1 , highest_vote_R1 = -1, highest_vote_R2 = -1, best_R2_i = -2, best_R2_j = -2;
+
+ for(j = 0; j < path_len ; j++){
+ if( this_masks[j] & CLUSTER_ALIGNMENT_DONOR_R1_MAPPED ){
+ if( highest_vote_R1 < vote_1 -> votes [ this_ii_path [j] ] [ this_jj_path [j] ] ){
+ best_R1_i = this_ii_path [j] ;
+ best_R1_j = this_jj_path [j] ;
+ highest_vote_R1 = vote_1 -> votes [ this_ii_path [j] ] [ this_jj_path [j] ] ;
+ }
+ }else{
+ if( highest_vote_R2 < vote_2 -> votes [ this_ii_path [j] ] [ this_jj_path [j] ] ){
+ best_R2_i = this_ii_path [j] ;
+ best_R2_j = this_jj_path [j] ;
+ highest_vote_R2 = vote_2 -> votes [ this_ii_path [j] ] [ this_jj_path [j] ] ;
+ }
+ }
+ //SUBREADprintf("MASK=%d\n", this_masks[j]);
+ }
+
+ //SUBREADprintf("IJ: R1=%d,%d R2=%d,%d MASK=%d\n", best_R1_i,best_R1_j,best_R2_i,best_R2_j);
+
+ simple_copy_vote_to_result( old_result_R1, vote_1, best_R1_i, best_R1_j , v1_all_subreads, vote_1 -> noninformative_subreads, this_score);
+ simple_copy_vote_to_result( old_result_R2, vote_2, best_R2_i, best_R2_j , v2_all_subreads, vote_2 -> noninformative_subreads, this_score);
+ break;
+ }
+ }
+ } else if( R1R2_mapped & ( CLUSTER_ALIGNMENT_DONOR_R1_MAPPED | CLUSTER_ALIGNMENT_DONOR_R2_MAPPED ) ) {
+ int is_R2_mapped = ( R1R2_mapped & CLUSTER_ALIGNMENT_DONOR_R2_MAPPED)?1:0;
+ for(i = 0; i < global_context -> config.multi_best_reads; i++) {
+ mapping_result_t * old_result_R = _global_retrieve_alignment_ptr(global_context, pair_number, is_R2_mapped, i);
+ short old_score_R = old_result_R -> selected_votes;
+
+
+
+ if( old_score_R < this_score ){
+
+ for(j = global_context -> config.multi_best_reads - 2; j>=i; j--){
+ mapping_result_t * shifted_result_R = _global_retrieve_alignment_ptr(global_context, pair_number, is_R2_mapped, j);
+ mapping_result_t * target_result_R = _global_retrieve_alignment_ptr(global_context, pair_number, is_R2_mapped, j+1);
+ memcpy( target_result_R, shifted_result_R , sizeof(mapping_result_t));
+
+ }
+
+ int best_R_i = -1, best_R_j = - 1 , highest_vote_R = -1;
+ gene_vote_t * this_vote = is_R2_mapped?vote_2:vote_1;
+
+ for(j = 0; j < path_len ; j++){
+ if( highest_vote_R < this_vote -> votes [ this_ii_path [j] ] [ this_jj_path [j] ] ){
+ best_R_i = this_ii_path [j] ;
+ best_R_j = this_jj_path [j] ;
+ highest_vote_R = this_vote -> votes [ this_ii_path [j] ] [ this_jj_path [j] ] ;
+ }
+ //SUBREADprintf("MASK=%d\n", this_masks[j]);
+ }
+
+ //SUBREADprintf("IJ: R1=%d,%d R2=%d,%d MASK=%d\n", best_R1_i,best_R1_j,best_R2_i,best_R2_j);
+
+ simple_copy_vote_to_result( old_result_R, this_vote, best_R_i, best_R_j , v1_all_subreads, this_vote -> noninformative_subreads, this_score);
+ break;
+ }
+
+
+ }
+ }
+ }
+
+ /*if(highest_score > 0){
+ if(this_score >0){
+ if(this_score > highest_score){
+ highest_score = this_score;
+ highest_occurance = 1;
+ memcpy(best_ii_path, this_ii_path, sizeof(int)*path_len);
+ memcpy(best_jj_path, this_jj_path, sizeof(int)*path_len);
+ best_path_len = path_len;
+ }else if(this_score == highest_score)
+ highest_occurance ++;
+ }
+ }*/
+ }
+
+
+ // call new junctions from the path
+ // then put the alignment into the best list.
+
+ }else{
+ return process_voting_junction_PE_topK(global_context, thread_context, pair_number, vote_1, vote_2, read_name_1, read_name_2, read_text_1, read_text_2, read_len_1, read_len_2, is_negative_strand, v1_all_subreads, v2_all_subreads);
+ }
+
+ free(cluster_buffer);
+ return 0;
+}
+
+
+int compare_cluster_elements (void * arr, int l, int r){
+ int * ii_array = ((void **)arr)[0];
+ int * jj_array = ((void **)arr)[1];
+ int * second_vote = ((void **)arr)[2];
+
+ if(second_vote[l] != second_vote[r])
+ return second_vote[l] - second_vote[r];
+
+ gene_vote_t * vote_1 = ((void **)arr)[3];
+ gene_vote_t * vote_2 = ((void **)arr)[4];
+
+
+ gene_vote_t * this_vote_L = second_vote[l]?vote_2:vote_1;
+ gene_vote_t * this_vote_R = second_vote[r]?vote_2:vote_1;
+
+ return this_vote_L->coverage_start[ii_array[l]][jj_array[l]] - this_vote_R -> coverage_start[ii_array[r]][jj_array[r]];
+}
+
+void exchange_cluster_elements (void * arr, int l, int r){
+ int * ii_array = ((void **)arr)[0];
+ int * jj_array = ((void **)arr)[1];
+ int * second_vote = ((void **)arr)[2];
+
+ int ti;
+ ti = ii_array[l];
+ ii_array[l] = ii_array[r];
+ ii_array[r]=ti;
+
+ ti = jj_array[l];
+ jj_array[l] = jj_array[r];
+ jj_array[r]=ti;
+
+ ti = second_vote[l];
+ second_vote[l] = second_vote[r];
+ second_vote[r]=ti;
+}
+
+#define NEW_EXTEND_SCAN_INTRON_LONGEST 6000
+#define NEW_EXTEND_SCAN_EXON_SHORTEST 14
+
+int find_path(global_context_t * global_context, thread_context_t * thread_context, int start_element_i, int target_element_i, int * ii_array, int * jj_array, int * is_second_vote_array, gene_vote_t * vote_1, gene_vote_t * vote_2, char * read_name_1, char * read_name_2, char * read_text_1, char * read_text_2, int read_len_1, int read_len_2, int is_negative_strand, int * this_mask , int * exon_last_base);
+int find_donor_receptor(global_context_t * global_context, thread_context_t * thread_context, char * rname, char * rtext, int rlen, int start_coverage, int end_coverage, unsigned int start_pos, unsigned int end_pos, int indels_in_start, int v1, int v2, int * misma_bases, int * matched_bases, int * is_negative_donor);
+int extend_uncovered_region_juncs(global_context_t * global_context, thread_context_t * thread_context, char * rname, char * rtext, int rlen, int scan_to_tail, unsigned int scan_chro_start, int scan_read_start, unsigned short expect_donor, int * mismatched_bases_after_start, int * first_exon_last_base, unsigned int * first_exon_first_base, int * ret_mismatched_bases, int * is_negative_donor){
+
+ if( scan_to_tail ) assert( scan_read_start < rlen - NEW_EXTEND_SCAN_EXON_SHORTEST );
+ else assert( scan_read_start >= NEW_EXTEND_SCAN_EXON_SHORTEST);
+
+ gene_value_index_t * value_index = thread_context?thread_context->current_value_index:global_context->current_value_index;
+ int x1, best_best_score = -1, best_best_occurance = 0;
+
+ unsigned long long matching_target = 0, rolling_bases = 0;
+
+ for(x1 = 0 ; x1 < 8 ; x1++){
+ int nch = scan_to_tail? rtext[ rlen - 2 - x1 ] : rtext[ 10 - x1 ] ;
+ matching_target = ( matching_target << 8 ) | nch;
+ }
+ if(0 && FIXLENstrcmp("V0112_0155:7:1101:13762:2349#ACTTGA", rname ) == 0 )
+ SUBREADprintf("TAG=%016llX\n",matching_target);
+
+ for(x1 = 0; x1 < NEW_EXTEND_SCAN_INTRON_LONGEST ; x1++){
+ int best_last_exon_base = -1, matched_in_the_uncovered_gap = -1, mismatched_bases = -1, extended_should_mismatch = -1;
+ unsigned int scan_cursor = scan_chro_start ;
+ if(scan_to_tail) scan_cursor+=x1;else scan_cursor-=x1;
+ unsigned long long nch = gvindex_get( value_index, scan_cursor );
+ if(scan_to_tail)
+ rolling_bases = (rolling_bases >> 8) | nch << 56;
+ else
+ rolling_bases = nch | ( rolling_bases << 8 );
+
+ //SUBREADprintf("MATCH:%016llX,%016llX\n", rolling_bases, matching_target);
+ if(rolling_bases == matching_target){
+ //SUBREADprintf("PNTT-M\n");
+ if(scan_to_tail) {
+ best_last_exon_base = find_donor_receptor(global_context, thread_context, rname, rtext, rlen, scan_read_start, rlen - 2 - 7, scan_chro_start, scan_cursor - (rlen - 2) , 0, 0,0, &mismatched_bases , &matched_in_the_uncovered_gap, is_negative_donor);
+ if(best_last_exon_base>0)
+ extended_should_mismatch = match_chro( rtext + best_last_exon_base , value_index, scan_chro_start + best_last_exon_base, rlen - best_last_exon_base, 0, global_context->config.space_type);
+ } else {
+ best_last_exon_base = find_donor_receptor(global_context, thread_context, rname, rtext, rlen, 10, scan_read_start, scan_cursor - 3 , scan_chro_start, 0, 0,0, &mismatched_bases , &matched_in_the_uncovered_gap,is_negative_donor);
+ if(best_last_exon_base>0)
+ extended_should_mismatch = match_chro( rtext, value_index, scan_chro_start, best_last_exon_base , 0, global_context->config.space_type);
+ }
+
+ }
+
+ if(best_last_exon_base >0 && extended_should_mismatch < ( scan_to_tail?( rlen - best_last_exon_base - 4 ):( best_last_exon_base - 4 )) && mismatched_bases < 2 ){
+ int this_score;
+ if(scan_to_tail) this_score = rlen - scan_read_start - mismatched_bases;
+ else this_score = scan_read_start - mismatched_bases;
+ if(best_best_score < this_score){
+ best_best_score = this_score;
+ (*mismatched_bases_after_start) = mismatched_bases;
+ (*first_exon_last_base) = best_last_exon_base;
+ (*first_exon_first_base) = scan_to_tail?( scan_cursor - (rlen - 2) ) : ( scan_cursor - 3 );
+ (*ret_mismatched_bases) = mismatched_bases;
+ best_best_occurance = 1;
+
+ }else if( best_best_score == this_score ) best_best_occurance ++;
+ }
+
+
+ if(0 && (!scan_to_tail) && best_last_exon_base >0 && extended_should_mismatch < best_last_exon_base - 4 && mismatched_bases < 2){
+ char out1pos[100], out2pos[100];
+ absoffset_to_posstr(global_context, scan_chro_start + best_last_exon_base, out1pos);
+ absoffset_to_posstr(global_context, scan_cursor - (rlen - 2) + best_last_exon_base, out2pos);
+ SUBREADprintf("LIMMISMA: %d < %d - 4\t\tfor %s\n" , extended_should_mismatch, best_last_exon_base ,rname);
+
+ SUBREADprintf("HEAD MATCH: %s - %s : MM=%d ; SPLIT=%d\t%s\n",out1pos, out2pos, mismatched_bases, best_last_exon_base, rname);
+
+ SUBREADprintf("R =%s\nS1=", rtext);
+ int x2;
+ for(x2 = 0; x2 < rlen; x2++){
+ if(x2 > best_last_exon_base + 16) SUBREADprintf(" ");
+ else{
+ int nch = gvindex_get( value_index, scan_cursor - 3 + x2 );
+ SUBREADprintf("%c", nch);
+ }
+ }
+ SUBREADprintf("\nS2=");
+
+ for(x2 = 0; x2 < rlen; x2++){
+ if(x2 > best_last_exon_base + 16 ) SUBREADprintf(" ");
+ else{
+ int nch = gvindex_get( value_index, scan_chro_start +x2);
+ SUBREADprintf("%c", nch);
+ }
+ }
+ SUBREADprintf("\n ");
+
+ for(x2 = 0; x2 < rlen; x2++){
+ if(x2 < best_last_exon_base ) SUBREADprintf(" ");
+ else if( x2 > best_last_exon_base + 1 ) SUBREADprintf(" ");
+ else SUBREADprintf("|");
+ }
+ SUBREADprintf("\n\n");
+ }
+ if(0 && scan_to_tail && best_last_exon_base >0 && extended_should_mismatch < rlen - best_last_exon_base - 4 && mismatched_bases < 2){
+
+ SUBREADprintf("LIMMISMA: %d < %d - 4\t\tfor %s\n" , extended_should_mismatch, (rlen - best_last_exon_base ),rname);
+ char out1pos[100], out2pos[100];
+ absoffset_to_posstr(global_context, scan_chro_start + best_last_exon_base, out1pos);
+ absoffset_to_posstr(global_context, scan_cursor - (rlen - 2) + best_last_exon_base, out2pos);
+ SUBREADprintf("TAIL MATCH: %s - %s : MM=%d ; SPLIT=%d\t%s\n",out1pos, out2pos, mismatched_bases, best_last_exon_base, rname);
+
+ SUBREADprintf("R =%s\nS1=", rtext);
+ int x2;
+ for(x2 = 0; x2 < rlen; x2++){
+ if(x2 < scan_read_start - 16) SUBREADprintf(" ");
+ else{
+ int nch = gvindex_get( value_index, x2 + scan_chro_start);
+ SUBREADprintf("%c", nch);
+ }
+ }
+ SUBREADprintf("\nS2=");
+
+ for(x2 = 0; x2 < rlen; x2++){
+ if(x2 < best_last_exon_base - 16 ) SUBREADprintf(" ");
+ else{
+ int nch = gvindex_get( value_index, scan_cursor - (rlen - 2)+x2);
+ SUBREADprintf("%c", nch);
+ }
+ }
+ SUBREADprintf("\n ");
+
+ for(x2 = 0; x2 < rlen; x2++){
+ if(x2 < best_last_exon_base ) SUBREADprintf(" ");
+ else if( x2 > best_last_exon_base + 1 ) SUBREADprintf(" ");
+ else SUBREADprintf("|");
+ }
+ SUBREADprintf("\n\n");
+ }
+ }
+ if(0&&best_best_occurance>0 && best_best_score>0)
+ SUBREADprintf("OCCR=%d : SCR=%d\n", best_best_occurance, best_best_score);
+ if(best_best_occurance == 1) return best_best_score;
+ return -1;
+}
+
+void simple_add_junction( global_context_t * global_context, thread_context_t * thread_context, unsigned int left_edge_wanted, unsigned int right_edge_wanted, int indel_at_junction, int is_negative_donors ){
+ char * chro_name_left, *chro_name_right;
+ int chro_pos_left,chro_pos_right;
+
+ locate_gene_position( left_edge_wanted , &global_context -> chromosome_table, &chro_name_left, &chro_pos_left);
+ locate_gene_position( right_edge_wanted , &global_context -> chromosome_table, &chro_name_right, &chro_pos_right);
+ if((! global_context->config.do_fusion_detection ) && chro_name_right!=chro_name_left) return;
+
+ //insert event
+ HashTable * event_table = NULL;
+ chromosome_event_t * event_space = NULL;
+ if(thread_context)
+ {
+ event_table = ((indel_thread_context_t *)thread_context -> module_thread_contexts[MODULE_INDEL_ID]) -> event_entry_table;
+ event_space = ((indel_thread_context_t *)thread_context -> module_thread_contexts[MODULE_INDEL_ID]) -> event_space_dynamic;
+ }
+ else
+ {
+ event_table = ((indel_context_t *)global_context -> module_contexts[MODULE_INDEL_ID]) -> event_entry_table;
+ event_space = ((indel_context_t *)global_context -> module_contexts[MODULE_INDEL_ID]) -> event_space_dynamic;
+ }
+ chromosome_event_t * search_return [MAX_EVENT_ENTRIES_PER_SITE];
+ chromosome_event_t * found = NULL;
+ int found_events = search_event(global_context, event_table, event_space, left_edge_wanted , EVENT_SEARCH_BY_SMALL_SIDE, CHRO_EVENT_TYPE_INDEL | CHRO_EVENT_TYPE_JUNCTION | CHRO_EVENT_TYPE_FUSION, search_return);
+
+ if(found_events)
+ {
+ int kx1;
+ for(kx1 = 0; kx1 < found_events ; kx1++)
+ {
+ if(search_return[kx1] -> event_large_side == right_edge_wanted)
+ {
+ found = search_return[kx1];
+ break;
+ }
+ }
+ }
+
+ if(found) found -> supporting_reads ++;
+ else
+ {
+ int event_no;
+
+
+ if(thread_context)
+ event_no = ((indel_thread_context_t *)thread_context -> module_thread_contexts[MODULE_INDEL_ID]) -> total_events ++;
+ else
+ event_no = ((indel_context_t *)global_context -> module_contexts[MODULE_INDEL_ID]) -> total_events ++;
+
+
+ chromosome_event_t * new_event = event_space+event_no;
+ memset(new_event,0,sizeof(chromosome_event_t));
+ new_event -> event_small_side = left_edge_wanted;
+ new_event -> event_large_side = right_edge_wanted;
+ new_event -> is_negative_strand= is_negative_donors;
+ new_event -> event_type = CHRO_EVENT_TYPE_JUNCTION;
+ new_event -> supporting_reads = 1;
+ new_event -> indel_length = 0;
+ new_event -> indel_at_junction = indel_at_junction;
+ new_event -> is_donor_found = 1;
+ new_event -> small_side_increasing_coordinate = 0;
+ new_event -> large_side_increasing_coordinate = 1;
+ put_new_event(event_table, new_event , event_no);
+ }
+}
+
+int align_cluster(global_context_t * global_context, thread_context_t * thread_context, struct cluster_element * this_cluster, char * read_name_1, char * read_name_2, char * read_text_1, char * read_text_2, int read_len_1, int read_len_2, int is_negative_strand, gene_vote_t * vote_1, gene_vote_t * vote_2, int * this_score, int * best_ii_path, int * best_jj_path, int * best_masks, int * best_path_length, int * R1R2_mapped){
+ //int cluster_x1;
+
+ //SUBREADprintf("\n === Cluster %s %s === \n", is_negative_strand?"NEG":"POS", read_name_1);
+ //unsigned int min_frag_start = 0xffffffff;
+
+ int ii_array[MAX_CLUSTER_ELEMENTS], jj_array[MAX_CLUSTER_ELEMENTS], is_second_vote_array[MAX_CLUSTER_ELEMENTS], dynamic_highest_mask[MAX_CLUSTER_ELEMENTS], x1;
+ void * sort_pointers [5];
+
+ for(x1 = 0 ; x1 < this_cluster->cluster_members; x1++){
+ ii_array[x1] = this_cluster -> i_list[x1];
+ jj_array[x1] = this_cluster -> j_list[x1];
+ is_second_vote_array[x1] = this_cluster -> from_second_read[x1];
+
+ }
+
+ sort_pointers[0] = ii_array;
+ sort_pointers[1] = jj_array;
+ sort_pointers[2] = is_second_vote_array;
+ sort_pointers[3] = vote_1;
+ sort_pointers[4] = vote_2;
+
+ basic_sort(sort_pointers, this_cluster->cluster_members, compare_cluster_elements, exchange_cluster_elements);
+
+ if(0)
+ for(x1 = 0 ; x1 < this_cluster->cluster_members; x1++){
+ gene_vote_t * this_vote = is_second_vote_array[x1]?vote_2:vote_1;
+ int ii = ii_array[x1];
+ int jj = jj_array[x1];
+
+ SUBREADprintf(" R%d %d - %d POS=%u VOTES=%d\n", 1+is_second_vote_array[x1], this_vote->coverage_start[ii][jj], this_vote->coverage_end[ii][jj], this_vote -> pos[ii][jj], this_vote->votes[ii][jj]);
+ }
+
+ int dynamic_highest_scores[MAX_CLUSTER_ELEMENTS], dynamic_last_exon[MAX_CLUSTER_ELEMENTS];
+ char dynamic_highest_path[MAX_CLUSTER_ELEMENTS];
+ memset(dynamic_highest_scores,0,sizeof(int)*MAX_CLUSTER_ELEMENTS);
+
+ int target_element_i;
+
+ for(target_element_i = 0; target_element_i < this_cluster->cluster_members; target_element_i++){
+ gene_vote_t * this_vote = is_second_vote_array[target_element_i]?vote_2:vote_1;
+ int ii = ii_array[target_element_i];
+ int jj = jj_array[target_element_i];
+ dynamic_highest_scores[target_element_i] = this_vote->coverage_end[ii][jj] - this_vote->coverage_start[ii][jj];
+ dynamic_highest_path[ target_element_i ] = -1;
+ }
+
+ int highest_score = -1;
+ int highest_target_end = -1;
+ for(target_element_i = 0; target_element_i < this_cluster->cluster_members; target_element_i++){
+ int start_element_i;
+ for(start_element_i = 0; start_element_i < this_cluster->cluster_members; start_element_i++){
+ if(target_element_i <= start_element_i) continue;
+ int this_mask = -1, breakpount_last_exon = -1;
+ int increasing_score = find_path(global_context, thread_context, start_element_i, target_element_i, ii_array, jj_array, is_second_vote_array, vote_1, vote_2, read_name_1, read_name_2, read_text_1, read_text_2, read_len_1, read_len_2, is_negative_strand, &this_mask, &breakpount_last_exon);
+ if(increasing_score >= 0 && increasing_score + dynamic_highest_scores[start_element_i] > dynamic_highest_scores[target_element_i]){
+ dynamic_highest_path[ target_element_i ] = start_element_i;
+ dynamic_highest_scores[target_element_i] = increasing_score + dynamic_highest_scores[start_element_i];
+ dynamic_highest_mask[ target_element_i ] = this_mask;
+ dynamic_last_exon[ target_element_i ] = breakpount_last_exon;
+ if( dynamic_highest_scores[target_element_i] > highest_score ){
+ highest_score = dynamic_highest_scores[target_element_i] ;
+ highest_target_end = target_element_i;
+ }
+ }
+ }
+ }
+
+
+ if(highest_target_end >=0 && highest_score > 160 - 159){
+ int is_on_path [MAX_CLUSTER_ELEMENTS];
+ memset(is_on_path,0,sizeof(int)*MAX_CLUSTER_ELEMENTS);
+
+ gene_vote_t * last_vote = is_second_vote_array[ highest_target_end ]?vote_2:vote_1;
+ int last_section_read_end = last_vote -> coverage_end[ ii_array [highest_target_end] ] [ jj_array [highest_target_end] ];
+ int this_rlen = is_second_vote_array[ highest_target_end ]?read_len_2 : read_len_1;
+ int this_votes = last_vote -> votes [ ii_array [highest_target_end] ] [ jj_array [highest_target_end] ];
+ int tail_first_exon_last_base_in_read=-1, tail_mismatched_bases=-1;
+ unsigned int tail_first_exon_first_base_on_chro, tail_mapped_section_pos;
+ int front_first_exon_last_base_in_read=-1, front_mismatched_bases=-1;
+ unsigned int front_first_exon_first_base_on_chro, front_mapped_section_pos ;
+ int front_score = 0, tail_score = 0, front_negative_donor = 0, tail_negative_donor = 0;
+
+ if(0 && last_section_read_end < this_rlen - NEW_EXTEND_SCAN_EXON_SHORTEST && this_votes > 4)
+ {
+ char * this_rname = is_second_vote_array[ highest_target_end ]?read_name_2:read_name_1;
+ char * this_rtext = is_second_vote_array[ highest_target_end ]?read_text_2:read_text_1;
+ int scan_to_tail = 1, mismatched_bases_after_start;
+ tail_mapped_section_pos = last_vote -> pos[ ii_array [highest_target_end] ] [ jj_array [highest_target_end] ] +
+ last_vote -> current_indel_cursor[ ii_array [highest_target_end] ] [ jj_array [highest_target_end] ] ;
+ if(0){
+ char out1pos[100];
+ absoffset_to_posstr(global_context, tail_mapped_section_pos, out1pos);
+ SUBREADprintf("RN=%s\nSTART=%u, READ_START=%d, READ_FIRTS_BASE_POS=%s\n", this_rname, tail_mapped_section_pos , last_section_read_end, out1pos);
+ }
+
+ tail_score = extend_uncovered_region_juncs(global_context, thread_context, this_rname, this_rtext , this_rlen, scan_to_tail, tail_mapped_section_pos, last_section_read_end , -1, & mismatched_bases_after_start, &tail_first_exon_last_base_in_read, &tail_first_exon_first_base_on_chro, &tail_mismatched_bases, &tail_negative_donor);
+ }
+ (*best_path_length) = 0;
+ if(highest_score>0) while(1){
+ best_ii_path[(*best_path_length)] = ii_array[highest_target_end];
+ best_jj_path[(*best_path_length)] = jj_array[highest_target_end];
+ best_masks[ (*best_path_length) ] = dynamic_highest_mask[highest_target_end];
+
+ if( dynamic_last_exon [ highest_target_end ] > 0 ) best_masks[ (*best_path_length) ] |= ( is_second_vote_array[ highest_target_end ]?CLUSTER_ALIGNMENT_DONOR_R2_MAPPED:CLUSTER_ALIGNMENT_DONOR_R1_MAPPED);
+
+ if( is_second_vote_array[ highest_target_end ] ) (*R1R2_mapped) |= CLUSTER_ALIGNMENT_DONOR_R2_MAPPED;
+ else (*R1R2_mapped) |= CLUSTER_ALIGNMENT_DONOR_R1_MAPPED;
+
+ (*best_path_length)++;
+
+ is_on_path[highest_target_end] = 1;
+ if( dynamic_highest_path[highest_target_end] == -1 ) break;
+ highest_target_end = dynamic_highest_path[highest_target_end];
+ }
+
+ gene_vote_t * first_vote = is_second_vote_array[ highest_target_end ]?vote_2:vote_1;
+ int first_section_read_start = first_vote -> coverage_start [ ii_array [highest_target_end] ] [ jj_array [highest_target_end] ] ;
+ this_votes = first_vote -> votes [ ii_array [highest_target_end] ] [ jj_array [highest_target_end] ];
+
+ if(0 && first_section_read_start > NEW_EXTEND_SCAN_EXON_SHORTEST && this_votes > 4){
+ char * this_rname = is_second_vote_array[ highest_target_end ]?read_name_2:read_name_1;
+ char * this_rtext = is_second_vote_array[ highest_target_end ]?read_text_2:read_text_1;
+ int scan_to_tail = 0, mismatched_bases_after_start;
+ front_mapped_section_pos = first_vote -> pos[ ii_array [highest_target_end] ] [ jj_array [highest_target_end] ];
+
+ front_score = extend_uncovered_region_juncs(global_context, thread_context, this_rname, this_rtext , this_rlen, scan_to_tail, front_mapped_section_pos, first_section_read_start , -1, & mismatched_bases_after_start, &front_first_exon_last_base_in_read, &front_first_exon_first_base_on_chro, &front_mismatched_bases, &front_negative_donor);
+
+ }
+
+ if(0 && front_score>0 && tail_score>0){
+
+ SUBREADprintf("\n>>> %s\n", read_name_1);
+
+ for(target_element_i = 0; target_element_i < this_cluster->cluster_members; target_element_i++){
+ SUBREADprintf("R%d\t", is_second_vote_array[target_element_i]+1);
+ }
+ SUBREADprintf("\n");
+
+ for(target_element_i = 0; target_element_i < this_cluster->cluster_members; target_element_i++){
+ gene_vote_t * this_vote = is_second_vote_array[target_element_i]?vote_2:vote_1;
+ int ii = ii_array[target_element_i];
+ int jj = jj_array[target_element_i];
+
+ SUBREADprintf("%d%c%d\t", this_vote->coverage_start[ii][jj], is_on_path[target_element_i]?'=':'-', this_vote->coverage_end[ii][jj]);
+ }
+ SUBREADprintf("\n");
+
+ for(target_element_i = 0; target_element_i < this_cluster->cluster_members; target_element_i++){
+ SUBREADprintf("%d\t", dynamic_highest_scores[target_element_i]);
+ }
+ SUBREADprintf("\n");
+ SUBREADprintf("Extra_scores = %d, %d\n", front_score, tail_score);
+ }
+
+ (*this_score) = highest_score + max(0, front_score) + max(0, tail_score);
+ int applied_score_cut=0;
+ if(((*R1R2_mapped) & CLUSTER_ALIGNMENT_DONOR_R1_MAPPED)&&( (*R1R2_mapped) & CLUSTER_ALIGNMENT_DONOR_R2_MAPPED ) )
+ applied_score_cut = read_len_2 + read_len_1 - 70;
+ else if((*R1R2_mapped) & CLUSTER_ALIGNMENT_DONOR_R1_MAPPED)
+ applied_score_cut = read_len_1 - 30;
+ else if((*R1R2_mapped) & CLUSTER_ALIGNMENT_DONOR_R1_MAPPED)
+ applied_score_cut = read_len_2 - 30;
+
+ if( (*this_score) >= applied_score_cut){
+ for( x1 = 0; x1 < MAX_CLUSTER_ELEMENTS; x1++){
+ if(!is_on_path[x1]) continue;
+
+ int x2, second_end = -1;
+ for(x2 = x1 + 1; x2 < MAX_CLUSTER_ELEMENTS; x2++){
+ if(is_on_path[x2]){
+ second_end = x2;
+ break;
+ }
+ }
+
+
+ if(second_end > 0){
+ if( dynamic_last_exon[second_end] >0 ){
+ gene_vote_t * this_vote = is_second_vote_array[ x1 ]?vote_2:vote_1;
+ unsigned int junction_small_side = this_vote -> pos[ ii_array[ x1 ]][ jj_array[ x1 ]] +
+ this_vote -> current_indel_cursor[ ii_array[ x1 ]][ jj_array[ x1 ]] + dynamic_last_exon[second_end];
+
+ unsigned int junction_large_side = this_vote -> pos[ ii_array[ second_end ]][ jj_array[ second_end ]] + dynamic_last_exon[second_end] + 1;
+
+ if(0){
+ char out1pos[100], out2pos[100];
+ absoffset_to_posstr(global_context, junction_small_side, out1pos);
+ absoffset_to_posstr(global_context, junction_large_side, out2pos);
+ SUBREADprintf("CLUSTER_JUNCTION %s %s\n%s\n%s\n\n", out1pos, out2pos, read_text_1, read_text_2);
+ }
+
+ //#warning "SUBREAD_151 ============= MAKE SURE: CHANGE '0' TO INSERTED BASES ================="
+ simple_add_junction(global_context, thread_context, junction_small_side, junction_large_side, 0, (dynamic_highest_mask[ second_end ] & CLUSTER_ALIGNMENT_DONOR_NEGATIVE_STRAND)?1:0);
+ }
+ }
+ }
+
+
+
+ if(0 && front_mismatched_bases <1 && front_score >14){
+ unsigned int junction_small_side = front_first_exon_first_base_on_chro + front_first_exon_last_base_in_read;
+ unsigned int junction_large_side = front_mapped_section_pos + front_first_exon_last_base_in_read + 1;
+
+ char out1pos[100], out2pos[100];
+ absoffset_to_posstr(global_context, junction_small_side+1, out1pos);
+ absoffset_to_posstr(global_context, junction_large_side+1, out2pos);
+ //SUBREADprintf("FMB=%d\tFS=%d\nPOS=%s , %s\n\n", front_mismatched_bases, front_score, out1pos, out2pos);
+
+ simple_add_junction(global_context, thread_context, junction_small_side, junction_large_side, 0, front_negative_donor);
+ }
+
+ if(0 && tail_mismatched_bases <1 && tail_score >14){
+ unsigned int junction_small_side = tail_mapped_section_pos + tail_first_exon_last_base_in_read;
+ unsigned int junction_large_side = tail_first_exon_first_base_on_chro + tail_first_exon_last_base_in_read;
+
+ char out1pos[100], out2pos[100];
+ absoffset_to_posstr(global_context, junction_small_side+1, out1pos);
+ absoffset_to_posstr(global_context, junction_large_side+1, out2pos);
+ //SUBREADprintf("BMB=%d\tBS=%d\nPOS=%s , %s\n\n", tail_mismatched_bases, tail_score, out1pos, out2pos);
+
+
+
+ simple_add_junction(global_context, thread_context, junction_small_side, junction_large_side, 0, tail_negative_donor);
+ }
+ }
+ }
+ return 0;
+}
+
+#define paired_donor_receptor_m2(s, c1, c2 ) ( s[0] == c1 && s[1] == c2 )
+
+int is_paired_donor_receptor( char * small_bases, char * large_bases ){
+
+ //SUBREADprintf("SITE1 = %c%c , SITE2 = %c%c\n", small_bases[0], small_bases[1], large_bases[0], large_bases[1]);
+ //
+ if ( paired_donor_receptor_m2( small_bases, 'G', 'T' ) &&
+ paired_donor_receptor_m2( large_bases, 'A', 'G' ) )
+ return 1;
+
+ if ( paired_donor_receptor_m2( small_bases, 'C', 'T' ) &&
+ paired_donor_receptor_m2( large_bases, 'A', 'C' ) )
+ return 2;
+
+ // http://www.ncbi.nlm.nih.gov/pmc/articles/PMC113136/
+ // the 99.24% of splice site pairs should be GT-AG,
+ // 0.69% GC-AG,
+ // 0.05% AT-AC
+ // and finally only 0.02% could consist of other types of non-canonical splice sites.
+
+ // non-canonical : GC-AG (+) or CT-GC (-)
+ if ( paired_donor_receptor_m2( small_bases, 'G', 'C' ) &&
+ paired_donor_receptor_m2( large_bases, 'A', 'G' ) )
+ return 3;
+
+ if ( paired_donor_receptor_m2( small_bases, 'C', 'T' ) &&
+ paired_donor_receptor_m2( large_bases, 'G', 'C' ) )
+ return 4;
+
+
+ // non-canonical : AT-AC (+) or GT-AT (-)
+ if ( paired_donor_receptor_m2( small_bases, 'A', 'T' ) &&
+ paired_donor_receptor_m2( large_bases, 'A', 'C' ) )
+ return 5;
+
+ if ( paired_donor_receptor_m2( small_bases, 'G', 'T' ) &&
+ paired_donor_receptor_m2( large_bases, 'A', 'T' ) )
+ return 6;
+
+
+ return 0;
+}
+
+int find_donor_receptor(global_context_t * global_context, thread_context_t * thread_context, char * rname, char * rtext, int rlen, int start_coverage, int end_coverage, unsigned int start_pos, unsigned int end_pos, int indels_in_start, int v1, int v2, int * misma_bases, int * matched_bases, int * is_negative_donor){
+
+ gene_value_index_t * value_index = thread_context?thread_context->current_value_index:global_context->current_value_index;
+ int search_in_read_start = start_coverage - 8, search_in_read_end = end_coverage + 8;
+ search_in_read_start = max(0, search_in_read_start);
+ search_in_read_end = min( rlen, search_in_read_end );
+ unsigned int search_in_chro_start = start_pos + indels_in_start + search_in_read_start;
+
+ char chro_bases_startside[ search_in_read_end - search_in_read_start ], chro_bases_endside[search_in_read_end - search_in_read_start];
+
+ int x1;
+
+ for(x1 = 0; x1 < search_in_read_end - search_in_read_start; x1++){
+ chro_bases_startside[x1] = gvindex_get( value_index, search_in_chro_start + x1 );
+ chro_bases_endside[x1] = gvindex_get( value_index , end_pos + search_in_read_start + x1);
+ }
+
+ int insertion_in_between_i, best_testing_score = 500 * 1000;
+ int best_insertion_in_between = -1, best_last_exon_base_in_start = -1;
+ int applied_insertion_limit = global_context->config.max_insertion_at_junctions;
+ for(insertion_in_between_i = 0; insertion_in_between_i <= applied_insertion_limit; insertion_in_between_i ++){
+ int start_site_match [ search_in_read_end - search_in_read_start ], end_site_match[ search_in_read_end - search_in_read_start ];
+ int start_last_exon_base, end_site_mismatches = 0, start_site_mismatches = 0;
+ for(start_last_exon_base = 0 ; start_last_exon_base < search_in_read_end - search_in_read_start ; start_last_exon_base++){
+ start_site_match[start_last_exon_base] = ( rtext[ search_in_read_start + start_last_exon_base ] == chro_bases_startside[start_last_exon_base] );
+ int end_site_x = ( rtext[ search_in_read_start + start_last_exon_base] == chro_bases_endside[start_last_exon_base] );
+ end_site_match[start_last_exon_base] = end_site_x;
+
+ if(start_last_exon_base >=insertion_in_between_i )
+ end_site_mismatches += !end_site_x;
+ }
+
+ for(start_last_exon_base = 0 ; start_last_exon_base < search_in_read_end - search_in_read_start - insertion_in_between_i ; start_last_exon_base++){
+ end_site_mismatches -= (! end_site_match[start_last_exon_base + insertion_in_between_i] );
+ start_site_mismatches += (! start_site_match[start_last_exon_base] );
+
+ if(start_last_exon_base >= 2 && start_last_exon_base < search_in_read_end - search_in_read_start -insertion_in_between_i -2){
+
+
+ // if(FIXLENstrcmp("V0112_0155:7:1101:12618:2466#ACTTGA", rname) == 0)
+ // SUBREADprintf("split=%d, ins=%d, MM=%d+%d \n", start_last_exon_base, insertion_in_between_i, start_site_mismatches, end_site_mismatches);
+
+
+ if( (end_site_mismatches + start_site_mismatches) * 500 + insertion_in_between_i < best_testing_score ){
+ int donor_paired_ret=is_paired_donor_receptor( chro_bases_startside + start_last_exon_base + 1, chro_bases_endside + insertion_in_between_i + start_last_exon_base - 1 );
+
+ if( donor_paired_ret ) {
+ best_insertion_in_between = insertion_in_between_i;
+ best_last_exon_base_in_start = start_last_exon_base;
+ best_testing_score = (end_site_mismatches + start_site_mismatches) * 500 + insertion_in_between_i;
+ (*misma_bases) = (end_site_mismatches + start_site_mismatches);
+ if(is_negative_donor) (*is_negative_donor) =(donor_paired_ret -1)%2;
+ (*matched_bases) = end_coverage - start_coverage - insertion_in_between_i - (end_site_mismatches + start_site_mismatches);
+ }
+ }
+
+ }
+ }
+ }
+
+
+ if(0 && FIXLENstrcmp("V0112_0155:7:1101:12618:2466", rname)==0)
+ {
+ chro_bases_startside[x1] = 0;
+ chro_bases_endside[x1] = 0;
+ char sp1s[200];
+ for(x1 =0; x1<200; x1++) sp1s[x1]=' ';
+ sp1s[search_in_read_start] =0;
+
+ char sp2s[200];
+ for(x1 =0; x1<200; x1++) sp2s[x1]=' ';
+ sp2s[search_in_read_start + best_insertion_in_between] =0;
+
+
+
+ char spE[200];
+ for(x1 =0; x1<200; x1++) spE[x1]=' ';
+ spE[search_in_read_start + best_last_exon_base_in_start] =0;
+
+
+ char spBB[200];
+ for(x1 =0; x1<200; x1++) spBB[x1]=' ';
+ spBB[ best_insertion_in_between] =0;
+
+
+ char out1pos[100];
+ absoffset_to_posstr(global_context, search_in_chro_start, out1pos);
+
+ if(1 || FIXLENstrcmp("chr14:105",out1pos)==0){
+ SUBREADprintf("POS=%s\t\tINS=%d\t\t%s\n", out1pos, best_insertion_in_between, rname);
+ SUBREADprintf("R= %s\nS1=%s%s\nS2=%s%s\n %s|%s|\n\n", rtext, sp1s, chro_bases_startside, sp1s, chro_bases_endside, spE, spBB);
+ }
+ }
+
+ if(best_last_exon_base_in_start>=0)
+ return best_last_exon_base_in_start + search_in_read_start ;
+ else return -1;
+}
+
+int find_path(global_context_t * global_context, thread_context_t * thread_context, int start_element_i, int target_element_i, int * ii_array, int * jj_array, int * is_second_vote_array, gene_vote_t * vote_1, gene_vote_t * vote_2, char * read_name_1, char * read_name_2, char * read_text_1, char * read_text_2, int read_len_1, int read_len_2, int is_negative_strand, int * this_mask , int * exon_last_base){
+ gene_vote_t * start_vote = is_second_vote_array[start_element_i]?vote_2:vote_1;
+ gene_vote_t * end_vote = is_second_vote_array[target_element_i]?vote_2:vote_1;
+
+ int start_coverage = start_vote->coverage_end[ ii_array[start_element_i] ][ jj_array[start_element_i] ];
+ int end_coverage = end_vote->coverage_start[ ii_array[target_element_i] ][ jj_array[target_element_i] ];
+ unsigned int start_pos = start_vote->pos[ ii_array[start_element_i] ][ jj_array[start_element_i] ];
+ unsigned int end_pos = end_vote->pos[ ii_array[target_element_i] ][ jj_array[target_element_i] ];
+ int ret = -1;
+
+ long long dist = start_pos;
+ dist -= end_pos;
+ (*this_mask)=0;
+ if( abs(dist)<50000 ) {
+ if(start_vote == end_vote){
+ if(start_coverage < end_coverage + 9){
+ char * this_read_name = is_second_vote_array[start_element_i]?read_name_2:read_name_1;
+ char * this_read_text = is_second_vote_array[start_element_i]?read_text_2:read_text_1;
+ int this_read_len = is_second_vote_array[start_element_i]?read_len_2:read_len_1, mismatched_bases = 0, matched_in_the_uncovered_gap = 0;
+ if(start_pos < end_pos){
+ int indels_in_start = start_vote -> current_indel_cursor [ ii_array[start_element_i]] [ jj_array[start_element_i] ] , donor_receptor_neg_strand = -1;
+ int best_last_base_in_start_exon = find_donor_receptor(global_context, thread_context, this_read_name, this_read_text, this_read_len, start_coverage, end_coverage, start_pos, end_pos, indels_in_start, start_vote -> votes[ ii_array[start_element_i]] [ jj_array[start_element_i] ], start_vote -> votes[ ii_array[target_element_i]] [ jj_array[target_element_i] ], &mismatched_bases , &matched_in_the_uncovered_gap, &donor_receptor_neg_strand);
+
+ if(best_last_base_in_start_exon > 0 && mismatched_bases<1){
+ ret = matched_in_the_uncovered_gap + end_vote->coverage_end[ ii_array[target_element_i] ][ jj_array[target_element_i] ] - end_coverage;
+ (*this_mask) = donor_receptor_neg_strand? CLUSTER_ALIGNMENT_DONOR_NEGATIVE_STRAND : 0 ;
+
+ if(0)SUBREADprintf("FROM %d-%d to %d-%d : INC=%d, UNCOV=%d/%d\n",
+ start_vote->coverage_start[ ii_array[start_element_i] ][ jj_array[start_element_i] ],
+ start_vote->coverage_end[ ii_array[start_element_i] ][ jj_array[start_element_i] ],
+ end_vote -> coverage_start[ ii_array[target_element_i] ][ jj_array[target_element_i] ],
+ end_vote -> coverage_end[ ii_array[target_element_i] ][ jj_array[target_element_i] ], ret,
+ matched_in_the_uncovered_gap , end_coverage - start_coverage);
+
+ // # of matched bases, from the end of the "start" section to the end of the end section.
+ *exon_last_base = best_last_base_in_start_exon;
+ }
+ }
+ }
+ }else{
+ ret = end_vote->coverage_end[ ii_array[target_element_i] ][ jj_array[target_element_i] ] - end_vote->coverage_start[ ii_array[target_element_i] ][ jj_array[target_element_i] ] ;
+ // if the two sections are on two reads, check the first base of the second read is after the first base of the first read.
+ }
+ }
+ return ret;
+}
+
int process_voting_junction_PE_topK(global_context_t * global_context, thread_context_t * thread_context, subread_read_number_t pair_number, gene_vote_t * vote_1, gene_vote_t * vote_2, char * read_name_1, char * read_name_2, char * read_text_1, char * read_text_2, int read_len_1, int read_len_2, int is_negative_strand, gene_vote_number_t v1_all_subreads, gene_vote_number_t v2_all_subreads)
{
vote_combination_t * comb_buffer = malloc(global_context -> config.max_vote_combinations * sizeof(vote_combination_t));
@@ -1434,7 +2222,9 @@ int process_voting_junction_PE_topK(global_context_t * global_context, thread_co
if(target_index < global_context -> config.max_vote_combinations){
int move_i;
+
for(move_i = min(used_comb_buffer, global_context -> config.max_vote_combinations - 1) ; move_i > target_index ; move_i --)
+ //checked: memory boundary
memcpy(comb_buffer + move_i, comb_buffer + move_i - 1 , sizeof(vote_combination_t) );
comb_buffer[target_index].r1_loc = vote_simple_1_buffer+i;
@@ -1519,8 +2309,10 @@ int process_voting_junction_PE_topK(global_context_t * global_context, thread_co
if(current_loc -> is_vote_t_item)
copy_vote_to_alignment_res(global_context, thread_context, current_alignment_tmp + (*current_r_cursor), current_junction_tmp ? current_junction_tmp + (*current_r_cursor) : NULL, current_vote, current_loc -> item_index_i, current_loc -> item_index_j, current_read_len, read_name_1, current_read_text, current_all_subreads , current_vote -> noninformative_subreads, pair_number, is_second_read, is_fully_covered);
else{
+ //checked: memory boundary
memcpy(current_alignment_tmp + (*current_r_cursor), _global_retrieve_alignment_ptr(global_context, pair_number, is_second_read, current_loc -> item_index_i), sizeof(mapping_result_t));
if(current_junction_tmp)
+ //checked: memory boundary
memcpy(current_junction_tmp + (*current_r_cursor), _global_retrieve_subjunc_ptr(global_context, pair_number, is_second_read, current_loc -> item_index_i), sizeof(subjunc_result_t));
}
(*current_r_cursor)++;
@@ -1573,8 +2365,10 @@ int process_voting_junction_PE_topK(global_context_t * global_context, thread_co
if(current_loc -> is_vote_t_item)
copy_vote_to_alignment_res(global_context, thread_context, current_alignment_tmp + (*current_r_cursor), current_junction_tmp ? current_junction_tmp + (*current_r_cursor): NULL, current_vote, current_loc -> item_index_i, current_loc -> item_index_j, current_read_len, read_name_1, current_read_text, current_all_subreads , current_vote -> noninformative_subreads, pair_number, is_second_read, is_fully_covered);
else{
+ //checked:boundary
memcpy(current_alignment_tmp + (*current_r_cursor), _global_retrieve_alignment_ptr(global_context, pair_number, is_second_read, current_loc -> item_index_i), sizeof(mapping_result_t));
if(current_junction_tmp)
+ //checked:boundary
memcpy(current_junction_tmp + (*current_r_cursor), _global_retrieve_subjunc_ptr(global_context, pair_number, is_second_read, current_loc -> item_index_i), sizeof(subjunc_result_t));
}
@@ -1593,6 +2387,14 @@ int process_voting_junction_PE_topK(global_context_t * global_context, thread_co
//SUBREADprintf("TOPK : CANDIDATES = %d , %d\n", alignment_res_r1_cursor, alignment_res_r2_cursor);
+ for(is_second_read = 0; is_second_read < 1 + global_context -> input_reads.is_paired_end_reads; is_second_read++){
+ int * current_r_cursor = is_second_read ? &alignment_res_r2_cursor:&alignment_res_r1_cursor;
+ if((*current_r_cursor) > global_context->config.multi_best_reads){
+ SUBREADprintf("ERROR: multi_best_locations excessed the boundary: %d > %d\n", (*current_r_cursor), global_context->config.multi_best_reads);
+ return -1;
+ }
+ }
+
for(is_second_read = 0; is_second_read < 1 + global_context -> input_reads.is_paired_end_reads; is_second_read++)
{
int * current_r_cursor = is_second_read ? &alignment_res_r2_cursor:&alignment_res_r1_cursor;
@@ -1605,6 +2407,7 @@ int process_voting_junction_PE_topK(global_context_t * global_context, thread_co
mapping_result_t * cur_res = _global_retrieve_alignment_ptr(global_context, pair_number, is_second_read, i);
if( i < (*current_r_cursor))
{
+ // checked:boundary
memcpy(cur_res, current_alignment_tmp + i, sizeof(mapping_result_t));
if(0 && FIXLENstrcmp("V0112_0155:7:1101:19612:13380", read_name_1)==0)
SUBREADprintf("COPIED READ_%d\t\t%llu [%d] , V=%d, MASK=%d, POS=%u, PTR=%p\n", is_second_read + 1, pair_number, *current_r_cursor, cur_res -> selected_votes, cur_res -> result_flags, current_alignment_tmp[i].selected_position, cur_res);
@@ -1615,6 +2418,7 @@ int process_voting_junction_PE_topK(global_context_t * global_context, thread_co
subjunc_result_t * cur_junc = _global_retrieve_subjunc_ptr(global_context, pair_number, is_second_read, i);
if(i < (*current_r_cursor))
{
+ // checked:boundary
memcpy(cur_junc, current_junction_tmp + i , sizeof(subjunc_result_t));
if(0 && FIXLENstrcmp("V0112_0155:7:1101:19612:13380", read_name_1)==0)
SUBREADprintf("COPIED SUBJUNC: MINOR=%u, MINORVOTES=%d\n", (current_junction_tmp + i) -> minor_position, (current_junction_tmp + i) -> minor_votes);
@@ -1712,10 +2516,11 @@ int is_funky_fragment(global_context_t * global_context, char * rname1, char * c
}
int process_voting_junction(global_context_t * global_context, thread_context_t * thread_context, subread_read_number_t pair_number, gene_vote_t * vote_1, gene_vote_t * vote_2, char * read_name_1, char * read_name_2, char * read_text_1, char * read_text_2, int read_len_1, int read_len_2, int is_negative_strand, gene_vote_number_t v1_all_subreads, gene_vote_number_t v2_all_subreads){
- //if(global_context -> input_reads.is_paired_end_reads || global_context -> config.do_breakpoint_detection)
- return process_voting_junction_PE_topK(global_context, thread_context, pair_number, vote_1, vote_2, read_name_1, read_name_2, read_text_1, read_text_2, read_len_1, read_len_2, is_negative_strand, v1_all_subreads, v2_all_subreads);
- //else
- // return process_voting_junction_SE(global_context, thread_context, pair_number, vote_1, read_name_1, read_text_1, read_len_1, is_negative_strand, v1_all_subreads);
+
+
+ //#warning "FOR TESTING CLUSTER_BASED JUNCTION DETECTION ONLY!!!"
+ //return process_voting_junction_PE_juncs(global_context, thread_context, pair_number, vote_1, vote_2, read_name_1, read_name_2, read_text_1, read_text_2, read_len_1, read_len_2, is_negative_strand, v1_all_subreads, v2_all_subreads);
+ return process_voting_junction_PE_topK(global_context, thread_context, pair_number, vote_1, vote_2, read_name_1, read_name_2, read_text_1, read_text_2, read_len_1, read_len_2, is_negative_strand, v1_all_subreads, v2_all_subreads);
}
@@ -1746,13 +2551,16 @@ unsigned int explain_read(global_context_t * global_context, thread_context_t *
explain_context.is_second_read = is_second_read ;
explain_context.best_read_id = best_read_id;
- unsigned int back_search_tail_position, front_search_start_position;
+ unsigned int back_search_tail_position,front_search_start_position;
unsigned short back_search_read_tail, front_search_read_start;
back_search_read_tail = min(explain_context.full_read_len , current_result -> confident_coverage_end );//- 5;
back_search_tail_position = current_result -> selected_position + back_search_read_tail + current_result -> indels_in_confident_coverage;
+ //if( back_search_read_tail > 102)
+ // SUBREADprintf("MAX back_search_read_tail : MIN %d , %d\n", explain_context.full_read_len , current_result -> confident_coverage_end);
+
explain_context.tmp_search_junctions[0].read_pos_end = back_search_read_tail;
explain_context.tmp_search_junctions[0].abs_offset_for_start = back_search_tail_position;
@@ -1772,8 +2580,8 @@ unsigned int explain_read(global_context_t * global_context, thread_context_t *
explain_context.best_is_pure_donor_found_explain = 0;
if(1) {
- front_search_read_start = back_search_read_tail - 8;
- front_search_start_position = back_search_tail_position - 8;
+ front_search_read_start = back_search_read_tail > 8? back_search_read_tail - 8:0;
+ front_search_start_position = back_search_tail_position>8?back_search_tail_position - 8:0;
} else {
//front_search_read_start = current_result -> confident_coverage_start + 5;
front_search_read_start = min(explain_context.full_read_len , current_result -> confident_coverage_end);
@@ -1847,9 +2655,15 @@ unsigned int explain_read(global_context_t * global_context, thread_context_t *
explain_context.tmp_search_junctions[0].abs_offset_for_start = front_search_start_position;
if(0 && FIXLENstrcmp("R000002689",explain_context.read_name ) == 0)
- SUBREADprintf("Enter F_SEARCH: start=%u read_pos=%d\n", front_search_start_position, front_search_read_start);
+ SUBREADprintf("Enter F_SEARCH: start=%u read_pos=%d REMAIN=%d\n", front_search_start_position, front_search_read_start, read_len - front_search_read_start );
- search_events_to_front(global_context, thread_context, &explain_context, read_text + front_search_read_start, qual_text + front_search_read_start, front_search_start_position,read_len - front_search_read_start , 0, 0, 1);
+
+ short search_remain = read_len - front_search_read_start;
+ //#warning "SUBREAD_151 REMOVE THE ASSERT! "
+ //if(search_remain >= 102)SUBREADprintf("FATAL: RLEN=%d, SEARCH=%d\n", read_len, front_search_read_start);
+ //assert( search_remain < 102 );
+
+ search_events_to_front(global_context, thread_context, &explain_context, read_text + front_search_read_start, qual_text + front_search_read_start, front_search_start_position, search_remain , 0, 0, 1);
if(0 && FIXLENstrcmp("R_chr901_932716_91M1D9M",explain_context.read_name ) == 0)
SUBREADprintf("F_SEARCH has found %d result sets\n", explain_context.all_front_alignments);
@@ -1919,11 +2733,9 @@ void debug_clipping(global_context_t * global_context, thread_context_t * threa
#define SOFT_CLIPPING_WINDOW_SIZE 5
#define SOFT_CLIPPING_MAX_ERROR 1
-#define find_soft_clipping_147 find_soft_clipping
-
// it returns the number of bases to be clipped off.
-int find_soft_clipping_147(global_context_t * global_context, thread_context_t * thread_context, gene_value_index_t * current_value_index, char * read_text, unsigned int mapped_pos, int test_len, int search_to_tail, int search_center)
+int find_soft_clipping(global_context_t * global_context, thread_context_t * thread_context, gene_value_index_t * current_value_index, char * read_text, unsigned int mapped_pos, int test_len, int search_to_tail, int search_center)
{
int base_in_window = 0;
int added_base_index = 0, removed_base_index = 0;
@@ -1938,7 +2750,7 @@ int find_soft_clipping_147(global_context_t * global_context, thread_context_t
else if(search_center >= test_len)
// SHOULD NOT HAPPEN!!!
search_start = test_len - 1;
- else search_start = search_center;
+ else search_start = search_center - 1;
delta = 1;
}else{
@@ -1947,7 +2759,7 @@ int find_soft_clipping_147(global_context_t * global_context, thread_context_t
search_start = 0;
else if(search_center >= test_len)
search_start = test_len - 1;
- else search_start = search_center;
+ else search_start = search_center + 1;
delta = -1;
}
@@ -1956,6 +2768,9 @@ int find_soft_clipping_147(global_context_t * global_context, thread_context_t
{
// add the new base
char reference_base = gvindex_get(current_value_index, added_base_index + mapped_pos);
+
+
+ //SUBREADprintf("CHMAT [%s] ref:read = %c:%c\n", search_to_tail?"T":"H", reference_base, read_text[added_base_index]);
int added_is_matched = (reference_base == read_text[added_base_index]);
matched_in_window += added_is_matched;
if(added_is_matched)
@@ -1993,53 +2808,10 @@ int find_soft_clipping_147(global_context_t * global_context, thread_context_t
else return search_start - 1;
}
}
-int find_soft_clipping_146(global_context_t * global_context, thread_context_t * thread_context, gene_value_index_t * current_value_index, char * read_text, unsigned int mapped_pos, int test_len, int search_to_tail, int search_center)
-{
-
- char window_matched[SOFT_CLIPPING_WINDOW_SIZE];
- int x0,x1,x2;
-
- memset(window_matched, 0 , SOFT_CLIPPING_WINDOW_SIZE);
-
- for(x0=0;x0 < test_len; x0++)
- {
-
- if(search_to_tail) x1 = test_len -1 -x0;
- else x1=x0;
- char ref_value = gvindex_get(current_value_index, mapped_pos + x1);
- int sum_matched=0;
- for(x2 = SOFT_CLIPPING_WINDOW_SIZE - 1; x2 > 0; x2--)
- {
- window_matched[x2] = window_matched[x2-1];
- sum_matched += window_matched[x2];
- }
- window_matched[0] = (ref_value == read_text[x1]);
- sum_matched += window_matched[0];
-
- /*
- for(x2 = 0; x2 < SOFT_CLIPPING_WINDOW_SIZE; x2++){
- printf("%d ", window_matched[x2]);
- }
- printf("\nMA=%d X0=%d\n", sum_matched, x0);
- */
-
- // find the first matched base, such that the matched bases >= SOFT_CLIPPING_WINDOW_SIZE - SOFT_CLIPPING_MAX_ERROR if this base is added into the window.
- if(window_matched[SOFT_CLIPPING_WINDOW_SIZE-1])
- {
- //SUBREADprintf("SOFTCLIP: %d > %d?\n", sum_matched, SOFT_CLIPPING_WINDOW_SIZE - SOFT_CLIPPING_MAX_ERROR);
- if(sum_matched >= SOFT_CLIPPING_WINDOW_SIZE - SOFT_CLIPPING_MAX_ERROR)
- {
- return max(0 , x0 - SOFT_CLIPPING_WINDOW_SIZE + 1);
- }
- }
-
- }
- return 0;
-}
// read_head_abs_offset is the first WANTED base in read.
// If the first section in read is reversed, read_head_abs_offset is the LAST WANTED bases in this section. (the abs offset of the first base in the section is actually larger than read_head_abs_offset)
-int final_CIGAR_quality(global_context_t * global_context, thread_context_t * thread_context, char * read_text, char * qual_text, int read_len, char * cigar_string, unsigned long read_head_abs_offset, int is_read_head_reversed, int * mismatched_bases, int covered_start, int covered_end, char * read_name, int * non_clipped_length, int *total_indel_length, int * matched_bases)
+int final_CIGAR_quality(global_context_t * global_context, thread_context_t * thread_context, char * read_text, char * qual_text, int read_len, char * cigar_string, unsigned long read_head_abs_offset, int is_read_head_reversed, int * mismatched_bases, int covered_start, int covered_end, char * read_name, int * non_clipped_length, int *total_indel_length, int * matched_bases, int * chromosomal_length)
{
int cigar_cursor = 0;
int read_cursor = 0;
@@ -2049,7 +2821,7 @@ int final_CIGAR_quality(global_context_t * global_context, thread_context_t * th
gene_value_index_t * current_value_index = thread_context?thread_context->current_value_index:global_context->current_value_index;
int current_reversed = is_read_head_reversed;
int all_mismatched = 0;
- int is_First_M = 1;
+ int is_First_M = 1, is_wrong_cigar = 0;
int head_soft_clipped = -1, tail_soft_clipped = -1;
unsigned int tmp_int = 0;
@@ -2062,6 +2834,8 @@ int final_CIGAR_quality(global_context_t * global_context, thread_context_t * th
if(isdigit(nch))
tmp_int = tmp_int*10+(nch-'0');
else{
+ if(tmp_int == 0)is_wrong_cigar = 1;
+ if(is_wrong_cigar) break;
if(nch == 'M' || nch == 'S')
{
char *qual_text_cur;
@@ -2080,6 +2854,7 @@ int final_CIGAR_quality(global_context_t * global_context, thread_context_t * th
int adj_coverage_start = covered_start - read_cursor;
char * debug_ptr = read_text;
+
if(current_reversed)
{
reversed_first_section_text = malloc(MAX_READ_LENGTH);
@@ -2110,6 +2885,7 @@ int final_CIGAR_quality(global_context_t * global_context, thread_context_t * th
if(current_reversed)
{
reversed_first_section_text = malloc(MAX_READ_LENGTH);
+ // checked: boundary
memcpy(reversed_first_section_text, read_text + read_cursor, tmp_int);
reverse_read(reversed_first_section_text, tmp_int, global_context->config.space_type);
debug_ptr = reversed_first_section_text;
@@ -2126,6 +2902,7 @@ int final_CIGAR_quality(global_context_t * global_context, thread_context_t * th
if(reversed_first_section_text)
free(reversed_first_section_text);
}
+
if(is_Last_M && is_First_M && tail_soft_clipped+head_soft_clipped >= tmp_int-1)
{
head_soft_clipped=0;
@@ -2181,6 +2958,11 @@ int final_CIGAR_quality(global_context_t * global_context, thread_context_t * th
if(nch == 'b') current_reversed = !current_reversed;
}
+ if(read_cursor>MAX_READ_LENGTH){
+ SUBREADprintf("ERROR: Cigar section longer than read length: %d >= %d, '%s'\n", tmp_int , MAX_READ_LENGTH, cigar_string);
+ is_wrong_cigar = 1;
+ }
+
tmp_int = 0;
}
}
@@ -2192,12 +2974,12 @@ int final_CIGAR_quality(global_context_t * global_context, thread_context_t * th
//#warning " ========== COMMENT THIS LINE !! ========="
//printf("QCR ALL MM=%d, RBLEN=%d, MAPPED_LEN=%d ; CIGAR=%s\n", all_mismatched, rebuilt_read_len , my_non_clipped_length, cigar_string);
- if(rebuilt_read_len != read_len || my_non_clipped_length < global_context->config.min_mapped_fraction){
+ if(is_wrong_cigar || rebuilt_read_len != read_len || my_non_clipped_length < global_context->config.min_mapped_fraction){
(*mismatched_bases)=99999;
all_matched_bases = 0;
sprintf(cigar_string, "%dM", read_len);
}
- else if(global_context -> config.show_soft_cliping && (head_soft_clipped>0 || tail_soft_clipped>0))
+ else if((head_soft_clipped>0 || tail_soft_clipped>0))
{
char new_cigar_tmp[120];
is_First_M=1;
@@ -2255,6 +3037,7 @@ int final_CIGAR_quality(global_context_t * global_context, thread_context_t * th
(*non_clipped_length) = my_non_clipped_length;
(*matched_bases) = my_non_clipped_length - all_mismatched - total_insertion_length;
+ (*chromosomal_length) = current_perfect_section_abs - read_head_abs_offset + total_insertion_length;
return max(0, (int)(all_matched_bases*60/my_non_clipped_length));
}
@@ -2280,9 +3063,14 @@ unsigned int finalise_explain_CIGAR(global_context_t * global_context, thread_co
// reverse the back_search result for every equally best alignment
//
for(back_i = 0; back_i < explain_context -> all_back_alignments; back_i++){
+ if( explain_context -> result_back_junction_numbers[back_i] > MAX_EVENTS_IN_READ ){
+ SUBREADprintf("ERROR: Too many cigar sections: %d > %d\n", explain_context -> result_back_junction_numbers[back_i] , MAX_EVENTS_IN_READ);
+ return 0;
+ }
for(xk1=0; xk1<explain_context -> result_back_junction_numbers[back_i]/2; xk1++)
{
perfect_section_in_read_t tmp_exp;
+ // checked: boundary
memcpy(&tmp_exp, &explain_context -> result_back_junctions[back_i][xk1], sizeof(perfect_section_in_read_t));
memcpy(&explain_context -> result_back_junctions[back_i][xk1], &explain_context -> result_back_junctions[back_i][explain_context -> result_back_junction_numbers[back_i] - xk1 - 1] , sizeof(perfect_section_in_read_t));
memcpy(&explain_context -> result_back_junctions[back_i][explain_context -> result_back_junction_numbers[back_i] - xk1 - 1] , &tmp_exp , sizeof(perfect_section_in_read_t));
@@ -2464,12 +3252,13 @@ unsigned int finalise_explain_CIGAR(global_context_t * global_context, thread_co
- int final_qual = 0, applied_mismatch = 0, non_clipped_length = 0, total_indel_length = 0, total_coverage_length = 0, final_MATCH = 0;
+ int final_qual = 0, applied_mismatch = 0, non_clipped_length = 0, total_indel_length = 0, total_coverage_length = 0, final_MATCH = 0, chromosomal_length = 0;
if(is_exonic_read_fraction_OK)
{
total_coverage_length = result -> confident_coverage_end - result -> confident_coverage_start;
- final_qual = final_CIGAR_quality(global_context, thread_context, explain_context -> full_read_text, explain_context -> full_qual_text, explain_context -> full_read_len , tmp_cigar, final_position, is_first_section_negative != ((result->result_flags & CORE_IS_NEGATIVE_STRAND)?1:0), &mismatch_bases, result -> confident_coverage_start, result -> confident_coverage_end, explain_context -> read_name, &non_clipped_length, &total_indel_length, & final_MATCH);
+ final_qual = final_CIGAR_quality(global_context, thread_context, explain_context -> full_read_text, explain_context -> full_qual_text, explain_context -> full_read_len , tmp_cigar, final_position, is_first_section_negative != ((result->result_flags & CORE_IS_NEGATIVE_STRAND)?1:0), &mismatch_bases, result -> confident_coverage_start, result -> confident_coverage_end, explain_context -> read_name, &non_clipped_length, &total_indel_length, & final_MATCH, & chromosomal_length);
+ //if(mismatch_bases<99) SUBREADprintf("CIGAR=%s, CHRLEN=%d\n", tmp_cigar, chromosomal_length);
applied_mismatch = is_junction_read? global_context->config.max_mismatch_junction_reads:global_context->config.max_mismatch_exonic_reads ;
@@ -2484,7 +3273,7 @@ unsigned int finalise_explain_CIGAR(global_context_t * global_context, thread_co
//if(explain_context -> pair_number == 999999)
// ACDB PVDB TTTS
- if(0 && FIXLENstrcmp("R001968841", explain_context -> read_name) ==0)
+ if(0 && FIXLENstrcmp("V0112_0155:7:1101:19274:15465", explain_context -> read_name) ==0)
SUBREADprintf("FINALQUAL %s : FINAL_POS=%u\tCIGAR=%s\tMM=%d > %d?\tVOTE=%d > %0.2f x %d ? MASK=%d\tQUAL=%d\tBRNO=%d\n\n", explain_context -> read_name, final_position , tmp_cigar, mismatch_bases, applied_mismatch, result -> selected_votes, global_context -> config.minimum_exonic_subread_fraction,result-> used_subreads_in_vote, result->result_flags, final_qual, explain_context -> best_read_id);
@@ -2496,6 +3285,7 @@ unsigned int finalise_explain_CIGAR(global_context_t * global_context, thread_co
realign_res -> realign_flags = result->result_flags;
realign_res -> first_base_is_jumpped = 0;
realign_res -> mapping_result = result;
+ realign_res -> chromosomal_length = chromosomal_length;
if(mismatch_bases > applied_mismatch ) realign_res -> realign_flags |= CORE_TOO_MANY_MISMATCHES;
else realign_res -> realign_flags &= ~CORE_TOO_MANY_MISMATCHES;
@@ -2962,50 +3752,11 @@ int donor_score(global_context_t * global_context, thread_context_t * thread_con
}
-#define NEW_EXTEND_SCAN_INTRON_LONGEST 5000
-#define NEW_EXTEND_SCAN_EXON_SHORTEST 12
-
-typedef struct {
- unsigned int small_exon_last_base;
- unsigned int large_exon_first_base;
- int canonical_donor_receptor_found;
-} newcore_extend_result_t;
-
-void newcore_extend_search_go(global_context_t * global_context, thread_context_t * thread_context, char * read_name, char * read_text, int search_to_tail, int candidate_last_base_in_exon_in_read, int candidate_last_base_in_exon_on_chro, newcore_extend_result_t * results, int * found_events) {
-
-}
-
-void newcore_extend_new_junctions( global_context_t * global_context, thread_context_t * thread_context, subread_read_number_t pair_number, char * read_name, char * read_text, char * qual_text, int read_len, int is_second_read, int best_read_id, mapping_result_t * result, subjunc_result_t * subjunc_result){
- int scan_to_tail;
- void * results;
- for(scan_to_tail = 0; scan_to_tail < 2 ; scan_to_tail++) {
- // (1) test if this read's worth scan to head and/or to tail
- int unexplained_head ;
- if(scan_to_tail) unexplained_head = read_len - result -> confident_coverage_end;
- else unexplained_head = result -> confident_coverage_start;
-
- if(unexplained_head < NEW_EXTEND_SCAN_EXON_SHORTEST) continue;
-
- // (2) scan to head or to tail
-
- unexplained_head += (scan_to_tail?-3:3);
- int candidate_last_base_in_exon_in_read = unexplained_head, found_events = 0;
- unsigned int candidate_last_base_in_exon_on_chro = result -> selected_position + unexplained_head;
-
- newcore_extend_search_go(global_context, thread_context, read_name, read_text, scan_to_tail, candidate_last_base_in_exon_in_read, candidate_last_base_in_exon_on_chro, results, &found_events);
- }
-}
-
-
void find_new_junctions(global_context_t * global_context, thread_context_t * thread_context, subread_read_number_t pair_number, char * read_name, char * read_text, char * qual_text, int read_len, int is_second_read, int best_read_id)
{
mapping_result_t * result =_global_retrieve_alignment_ptr(global_context, pair_number, is_second_read, best_read_id);
subjunc_result_t * subjunc_result =_global_retrieve_subjunc_ptr(global_context, pair_number, is_second_read, best_read_id);
-
- if(0)
- newcore_extend_new_junctions(global_context, thread_context, pair_number, read_name, read_text, qual_text, read_len, is_second_read, best_read_id, result, subjunc_result);
-
if(read_len > EXON_LONG_READ_LENGTH)
{
assert(result -> selected_position <= 0xffff0000);
@@ -3024,7 +3775,7 @@ void find_new_junctions(global_context_t * global_context, thread_context_t * th
{
- if(0 && FIXLENstrcmp("R006856515", read_name) == 0 )
+ if(0 && FIXLENstrcmp("R003738400", read_name) == 0 )
{
char posout[100];
int xk1;
@@ -3173,7 +3924,7 @@ void find_new_junctions(global_context_t * global_context, thread_context_t * th
// note that selected_real_split_point is the first UNWANTED base after left half.
//if(abs(left_edge_wanted-27286396) < 250 || abs(right_edge_wanted - 27286396)<250)
- if(0 && FIXLENstrcmp("V0112_0155:7:1101:19612:13380", read_name) == 0)
+ if(0 && FIXLENstrcmp("R003738400", read_name) == 0)
{
char leftpos[100], rightpos[100];
absoffset_to_posstr(global_context, left_edge_wanted, leftpos);
@@ -3374,8 +4125,9 @@ int write_fusion_final_results(global_context_t * global_context)
chro_pos_left++;
all_junctions ++;
- fprintf(ofp, "%s\t%u\t%s\t%u\t%s\t%d\n", chro_name_left, chro_pos_left, chro_name_right, chro_pos_right+1, event_body -> is_strand_jumped?"No":"Yes", event_body -> final_counted_reads);
- //fprintf(ofp, "%s\t%u\t%s\t%u\t%s\t%d\t%s\t%s\n", chro_name_left, chro_pos_left, chro_name_right, chro_pos_right+1, event_body -> is_strand_jumped?"No":"Yes", event_body -> final_counted_reads, event_body -> small_side_increasing_coordinate?"Yes":"No", event_body -> large_side_increasing_coordinate?"Yes":"No");
+ //#warning "SUBREAD_151 ================ COMMENT the 'UNPAIRED' line and UNCOMMENT the next line ======================"
+ //fprintf(ofp, "UNPAIRED\t%s\t%u\t%s\t%u\t%s\t%d\n", chro_name_left, chro_pos_left, chro_name_right, chro_pos_right, event_body -> is_strand_jumped?"No":"Yes", event_body -> final_counted_reads);
+ fprintf(ofp, "%s\t%u\t%s\t%u\t%s\t%d\t%s\t%s\n", chro_name_left, chro_pos_left, chro_name_right, chro_pos_right+1, event_body -> is_strand_jumped?"No":"Yes", event_body -> final_counted_reads, event_body -> small_side_increasing_coordinate?"Yes":"No", event_body -> large_side_increasing_coordinate?"Yes":"No");
}
global_context -> all_fusions = all_junctions;
diff --git a/src/core.c b/src/core.c
index 4a2c4f3..7423c0e 100644
--- a/src/core.c
+++ b/src/core.c
@@ -1124,6 +1124,60 @@ unsigned int move_to_read_head(unsigned int tailpos, char * cigar){
return tailpos;
}
+
+// This function returns 1 if the cut was added.
+// It returns 0 if the head or tail cut is not able to be added (e.g., the cigar ends up like "50M4I10S" or "10S30N90M")
+int add_head_tail_cut_softclipping(char * cigar, int rlen, int head_cut, int tail_cut){
+
+ char cigar_added [CORE_MAX_CIGAR_STR_LEN];
+ int cigar_cursor = 0, read_cursor = 0, next_read_cursor = 0;
+ int tmpi = 0, nch, has_M = 0;
+
+ cigar_added[0]=0;
+
+ while(1){
+ nch = cigar[cigar_cursor++];
+ if(nch == 0) break;
+
+ if(isdigit(nch)){
+ tmpi = 10 * tmpi + (nch - '0');
+ }else{
+ if('M' == nch || 'S' == nch || 'I' == nch)
+ next_read_cursor = read_cursor + tmpi;
+
+ int head_S = 0, tail_S =0, remainder_tmpi = tmpi, is_skip = 0;
+
+ if(next_read_cursor <= head_cut) is_skip = 1;
+ if(read_cursor >= rlen - tail_cut) is_skip = 1;
+ if(!is_skip){
+ if(nch != 'S'){
+ if(read_cursor <= head_cut)head_S = head_cut;
+ if(next_read_cursor >= rlen - tail_cut) tail_S = tail_cut;
+ remainder_tmpi = tmpi;
+ if(head_S) remainder_tmpi-= (head_S - read_cursor);
+ if(tail_S) remainder_tmpi-= next_read_cursor - (rlen - tail_S);
+ }
+ //SUBREADprintf("OPTR=%c, THREE-LENS=%d,%d,%d\n", nch, head_S, remainder_tmpi, tail_S);
+ if(head_S >0 || tail_S > 0) if(nch !='M') return 0;
+
+ if(head_S > 0) sprintf(cigar_added + strlen(cigar_added), "%dS", head_S );
+ if(remainder_tmpi > 0){
+ sprintf(cigar_added + strlen(cigar_added), "%d%c", remainder_tmpi , nch );
+ if( nch == 'M' ) has_M = 1;
+ }
+ if(tail_S > 0) sprintf(cigar_added + strlen(cigar_added), "%dS", tail_S );
+ }
+
+ read_cursor = next_read_cursor;
+ tmpi = 0;
+ }
+ }
+ //SUBREADprintf("RPCIGAR: %s => %s\n", cigar, cigar_added);
+ strcpy(cigar, cigar_added);
+ return has_M;
+}
+
+
int convert_read_to_tmp(global_context_t * global_context , subread_output_context_t * output_context, int read_number, int is_second_read, int read_len, char * read_text, char * qual_text, realignment_result_t * current_result, subread_output_tmp_t * r, char * read_name)
{
@@ -1180,7 +1234,7 @@ int convert_read_to_tmp(global_context_t * global_context , subread_output_conte
//if(( r->out_strands[xk1] == '-' ) != (r->out_strands[0]=='-' )) vitual_head_pos = move_to_read_head(vitual_head_pos, r->out_cigars[xk1]);
- if(0==locate_gene_position_max(vitual_head_pos ,& global_context -> chromosome_table, & chimaric_chr, & chimeric_pos, 0+r->out_lens[xk1]))
+ if(0==locate_gene_position_max(vitual_head_pos ,& global_context -> chromosome_table, & chimaric_chr, & chimeric_pos, NULL, NULL, 0+r->out_lens[xk1]))
{
int soft_clipping_movement = 0;
soft_clipping_movement = get_soft_clipping_length(r->out_cigars[xk1]);
@@ -1199,14 +1253,35 @@ int convert_read_to_tmp(global_context_t * global_context , subread_output_conte
if(is_r_OK)
{
- if(locate_gene_position_max(r->linear_position,& global_context -> chromosome_table, &r-> chro , &r -> offset, read_len))
- {
+ int head_cut = 0 , tail_cut = 0;
+
+ if(0 && FIXLENstrcmp("V0112_0155:7:1302:9507:32993", read_name)==0){
+ char posout1[100];
+ absoffset_to_posstr(global_context, r->linear_position, posout1);
+ SUBREADprintf("PERR : CIGAR=%s, READLEN=%d, POS=%s\n", r->cigar , read_len, posout1);
+ }
+
+ if(locate_gene_position_max(r->linear_position,& global_context -> chromosome_table, &r-> chro , &r -> offset, &head_cut, &tail_cut, global_context->config.do_fusion_detection?read_len:current_result->chromosomal_length)) {
is_r_OK = 0;
+ } else {
+
+
+ if(0 && FIXLENstrcmp("V0112_0155:7:1302:9507:32993", read_name)==0){
+ char posout1[100];
+ absoffset_to_posstr(global_context, r->linear_position, posout1);
+ SUBREADprintf("CUTT : CIGAR=%s, READLEN=%d, CATS=%d %d\n", r->cigar , read_len, head_cut, tail_cut);
}
- else
- {
- r -> offset++;
- assert(r-> chro);
+
+ int is_added_OK = 1;
+ if(head_cut!=0 || tail_cut!=0)
+ is_added_OK = add_head_tail_cut_softclipping(r->cigar , read_len, head_cut, tail_cut);
+
+ if(is_added_OK){
+ r -> offset++;
+ assert(r-> chro);
+ }else{
+ is_r_OK = 0;
+ }
}
if(global_context -> config.do_breakpoint_detection && !(current_result -> realign_flags & CORE_NOTFOUND_DONORS))
@@ -1718,21 +1793,23 @@ void write_single_fragment(global_context_t * global_context, thread_context_t *
if(is_R1_OK && is_R2_OK && rec1->chro == rec2->chro)tlen = calc_tlen(global_context , rec1, rec2, read_len_1, read_len_2);
+ //#warning "SUBREAD_151 ============= REMOVE '+1' when not doing structure variance detection for the proposal =========="
if(0)
if(global_context -> config.do_fusion_detection && rec1 && rec2 && current_location == 0 && rec1 -> chimeric_sections == 1 && rec2 -> chimeric_sections == 1){
// when current_location == 0, both r1 and r2 were on the same strand.
int is_funky = is_funky_fragment(global_context, read_name_1, rec1 -> chro, rec1 -> offset, read_len_1, rec1 -> strand, rec1 -> out_cigars[0], read_text_1, read_name_2, rec2 -> chro, rec2 -> offset, read_len_2, rec2 -> strand, rec2 -> out_cigars[0] , read_text_2, tlen);
- //if (is_funky)
- // SUBREADprintf("RNAME %s is %d, CHRO = '%s' %p '%s' %p, POS=%u, %u\n", read_name_1, is_funky, rec1 -> chro,rec1 -> chro,rec2 -> chro,rec2 -> chro, rec1 -> offset, rec2 -> offset);
+ if (0 &&is_funky)
+ SUBREADprintf("RNAME %s is %d, CHRO = '%s' %p '%s' %p, POS=%u, %u\n", read_name_1, is_funky, rec1 -> chro,rec1 -> chro,rec2 -> chro,rec2 -> chro, rec1 -> offset, rec2 -> offset);
+
if(is_funky & FUNKY_FRAGMENT_A){
fraglist_append(&global_context -> funky_list_A, pair_number);
}
- if(is_funky & FUNKY_FRAGMENT_BC){
+ if(1)if(is_funky & FUNKY_FRAGMENT_BC){
//#warning "LOGIC WRONG: R1 AND R2 SHOULD BE DECIDED BY THEIR MAPPING POSITIONS"
bktable_append(&global_context -> funky_table_BC, rec1 -> chro, rec1 -> offset + rec1 -> soft_clipping_movements, NULL + (2*pair_number));
bktable_append(&global_context -> funky_table_BC, rec2 -> chro, rec2 -> offset + rec2 -> soft_clipping_movements, NULL + (2*pair_number+1));
}
- if(is_funky & FUNKY_FRAGMENT_DE){
+ if(1)if(is_funky & FUNKY_FRAGMENT_DE){
fraglist_append(&global_context -> funky_list_DE, pair_number);
bktable_append(&global_context -> funky_table_DE, rec1 -> chro, rec1 -> offset + rec1 -> soft_clipping_movements, NULL + (2*pair_number + (rec1 -> offset > rec2 -> offset ? 1:0)));
bktable_append(&global_context -> funky_table_DE, rec2 -> chro, rec2 -> offset + rec2 -> soft_clipping_movements, NULL + (2*pair_number + (rec1 -> offset < rec2 -> offset ? 1:0)));
@@ -2140,7 +2217,8 @@ int do_iteration_one(global_context_t * global_context, thread_context_t * threa
sqr_read_number = 0;
}
}
- bigtable_release_result(global_context, thread_context, current_read_number, 1);
+
+ //bigtable_release_result(global_context, thread_context, current_read_number, 1);
}
@@ -2759,7 +2837,7 @@ int do_iteration_two(global_context_t * global_context, thread_context_t * threa
sqr_read_number=0;
}
}
- bigtable_release_result(global_context, thread_context, current_read_number, 1);
+ //bigtable_release_result(global_context, thread_context, current_read_number, 1);
}
free(final_realignments);
@@ -2860,7 +2938,7 @@ int do_voting(global_context_t * global_context, thread_context_t * thread_conte
int is_reversed, applied_subreads = 0, v1_all_subreads=0, v2_all_subreads=0;
ret = fetch_next_read_pair(global_context, thread_context, ginp1, ginp2, &read_len_1, &read_len_2, read_name_1, read_name_2, read_text_1, read_text_2, qual_text_1, qual_text_2,1, ¤t_read_number);
- //SUBREADprintf("DO_VOTE:%llu\n", current_read_number);
+ //SUBREADprintf("DO_VOTE:%llu BY THREAD %d\n", current_read_number, thread_context -> thread_id);
if(current_read_number < 0) break;
//SUBREADprintf("RL=%d,%d\n", read_len_1, read_len_2);
@@ -2987,8 +3065,7 @@ int do_voting(global_context_t * global_context, thread_context_t * thread_conte
{
//SUBREADprintf("P%d %llu %s\n", is_reversed, current_read_number, read_name_1);
- if(0 && FIXLENstrcmp("R001968841", read_name_1) ==0 )
- {
+ if(0 && FIXLENstrcmp("R003577537", read_name_1) ==0 ) {
SUBREADprintf(">>>%llu<<<\n%s [%d] %s\n%s [%d] %s\n", current_read_number, read_name_1, read_len_1, read_text_1, read_name_2, read_len_2, read_text_2);
SUBREADprintf(" ======= PAIR %s = %llu ; NON_INFORMATIVE = %d, %d =======\n", read_name_1, current_read_number, vote_1 -> noninformative_subreads, vote_2 -> noninformative_subreads);
print_votes(vote_1, global_context -> config.index_prefix);
@@ -2999,7 +3076,7 @@ int do_voting(global_context_t * global_context, thread_context_t * thread_conte
/*
if(407229 == current_read_number){
fprintf(stderr, "TABLE_ITEMS=%llu\n", global_context->current_index->current_items);
- print_votes(vote_1, global_context -> config.index_prefix);
+ //print_votes(vote_1, global_context -> config.index_prefix);
//print_votes(vote_2, global_context -> config.index_prefix);
}*/
//if(global_context -> input_reads.is_paired_end_reads) finalise_vote(vote_2);
@@ -3046,8 +3123,6 @@ int do_voting(global_context_t * global_context, thread_context_t * thread_conte
SUBREADprintf("FIN R %s : V=%d , %d",read_name_1, current_result_1 -> selected_votes, current_result_2 -> selected_votes);
}
- if(is_reversed==1)
- bigtable_release_result(global_context, thread_context, current_read_number, 1);
}
@@ -3107,6 +3182,7 @@ int do_voting(global_context_t * global_context, thread_context_t * thread_conte
}
+ //bigtable_release_result(global_context, thread_context, current_read_number, 1);
sqr_read_number++;
}
@@ -3134,6 +3210,7 @@ void * run_in_thread(void * pthread_param)
if(thread_init_lock)
subread_lock_release(thread_init_lock);
+
switch(task)
{
case STEP_VOTING:
@@ -3165,12 +3242,14 @@ int run_maybe_threads(global_context_t *global_context, int task)
void * thr_parameters [5];
int ret_value =0;
+ /*
if(task==STEP_VOTING)
print_in_box(80,0,0, "Map %s...", global_context->input_reads.is_paired_end_reads?"fragments":"reads");
else if(task == STEP_ITERATION_ONE)
print_in_box(80,0,0, "Detect indels%s...", global_context->config.do_breakpoint_detection?" and junctions":"");
else if(task == STEP_ITERATION_TWO)
print_in_box(80,0,0, "Finish the %'llu %s...", global_context -> processed_reads_in_chunk, global_context->input_reads.is_paired_end_reads?"fragments":"reads");
+ */
if(global_context->config.all_threads<2)
{
@@ -3480,55 +3559,51 @@ int print_configuration(global_context_t * context)
{
if(context->config.do_fusion_detection)
{
- print_in_box(80, 0, 0, " Function : Read alignment + Junction/Fusion detection%s", context->config.experiment_type == CORE_EXPERIMENT_DNASEQ?" (DNA-Seq)":" (RNA-Seq)");
+ print_in_box(80, 0, 0, "Function : Read alignment + Junction/Fusion detection%s", context->config.experiment_type == CORE_EXPERIMENT_DNASEQ?" (DNA-Seq)":" (RNA-Seq)");
}
else
- print_in_box(80, 0, 0, " Function : Read alignment + Junction detection (%s)", context->config.experiment_type == CORE_EXPERIMENT_DNASEQ?"DNA-Seq":"RNA-Seq");
+ print_in_box(80, 0, 0, "Function : Read alignment + Junction detection (%s)", context->config.experiment_type == CORE_EXPERIMENT_DNASEQ?"DNA-Seq":"RNA-Seq");
}
else
- print_in_box(80, 0, 0, " Function : Read alignment%s", context->config.experiment_type == CORE_EXPERIMENT_DNASEQ?" (DNA-Seq)":" (RNA-Seq)");
- print_in_box(80, 0, 0, " Threads : %d", context->config.all_threads);
+ print_in_box(80, 0, 0, "Function : Read alignment%s", context->config.experiment_type == CORE_EXPERIMENT_DNASEQ?" (DNA-Seq)":" (RNA-Seq)");
if( context->config.second_read_file[0])
{
- print_in_box(80, 0, 0, " Input file 1 : %s", context->config.first_read_file);
- print_in_box(80, 0, 0, " Input file 2 : %s", context->config.second_read_file);
+ print_in_box(80, 0, 0, "Input file 1 : %s", context->config.first_read_file);
+ print_in_box(80, 0, 0, "Input file 2 : %s", context->config.second_read_file);
}
else
- print_in_box(80, 0, 0, " Input file : %s%s", context->config.first_read_file, context->config.is_SAM_file_input?(context->config.is_BAM_input?" (BAM)":" (SAM)"):"");
+ print_in_box(80, 0, 0, "Input file : %s%s", context->config.first_read_file, context->config.is_SAM_file_input?(context->config.is_BAM_input?" (BAM)":" (SAM)"):"");
if(context->config.output_prefix [0])
- print_in_box(80, 0, 0, " Output file : %s (%s)", context->config.output_prefix, context->config.is_BAM_output?"BAM":"SAM");
+ print_in_box(80, 0, 0, "Output file : %s (%s)", context->config.output_prefix, context->config.is_BAM_output?"BAM":"SAM");
else
- print_in_box(80, 0, 0, " Output method : STDOUT (%s)" , context->config.is_BAM_output?"BAM":"SAM");
+ print_in_box(80, 0, 0, "Output method : STDOUT (%s)" , context->config.is_BAM_output?"BAM":"SAM");
- print_in_box(80, 0, 0, " Index name : %s", context->config.index_prefix);
- print_in_box(80, 0, 0, " Phred offset : %d", (context->config.phred_score_format == FASTQ_PHRED33)?33:64);
- //print_in_box(80, 0, 0, " Space type : %s", (context->config.space_type == GENE_SPACE_COLOR)?"color-space":"base-space");
- print_in_box(80, 0, 1, "");
+ print_in_box(80, 0, 0, "Index name : %s", context->config.index_prefix);
+ print_in_box(80, 0, 0, "");
+ print_in_box(80, 0, 0, " Threads : %d", context->config.all_threads);
+ print_in_box(80, 0, 0, " Phred offset : %d", (context->config.phred_score_format == FASTQ_PHRED33)?33:64);
if( context->config.second_read_file[0])
{
- print_in_box(80, 0, 0, " All subreads : %d", context->config.total_subreads);
- print_in_box(80, 0, 0, " Min read1 votes : %d", context->config.minimum_subread_for_first_read);
- print_in_box(80, 0, 0, " Min read2 votes : %d", context->config.minimum_subread_for_second_read);
- print_in_box(80, 0, 0, " Max fragment size : %d", context->config.maximum_pair_distance);
- print_in_box(80, 0, 0, " Min fragment size : %d", context->config.minimum_pair_distance);
- print_in_box(80, 0, 1, "");
+ print_in_box(80, 0, 0, " # of extracted subreads : %d", context->config.total_subreads);
+ print_in_box(80, 0, 0, " Min read1 vote : %d", context->config.minimum_subread_for_first_read);
+ print_in_box(80, 0, 0, " Min read2 vote : %d", context->config.minimum_subread_for_second_read);
+ print_in_box(80, 0, 0, " Max fragment size : %d", context->config.maximum_pair_distance);
+ print_in_box(80, 0, 0, " Min fragment size : %d", context->config.minimum_pair_distance);
}
else
- print_in_box(80, 0, 0, " Min votes : %d / %d", context->config.minimum_subread_for_first_read, context->config.total_subreads);
+ print_in_box(80, 0, 0, " Min votes : %d / %d", context->config.minimum_subread_for_first_read, context->config.total_subreads);
- print_in_box(80, 0, 0, " Allowed mismatch : %d bases", context->config.max_mismatch_exonic_reads);
- print_in_box(80, 0, 0, " Max indels : %d", context->config.max_indel_length);
- print_in_box(80, 0, 0, " # of Best mapping : %d", context->config.multi_best_reads);
- print_in_box(80, 0, 0, " Unique mapping : %s", context->config.report_multi_mapping_reads?"no":"yes");
- print_in_box(80, 0, 0, " Hamming distance : %s", context->config.use_hamming_distance_break_ties?"yes":"no");
- print_in_box(80, 0, 0, " Quality scores : %s", context->config.use_quality_score_break_ties?"yes":"no");
+ print_in_box(80, 0, 0, " Maximum allowed mismatches : %d", context->config.max_mismatch_exonic_reads);
+ print_in_box(80, 0, 0, " Maximum allowed indel bases : %d", context->config.max_indel_length);
+ print_in_box(80, 0, 0, "# of best alignments reported : %d", context->config.multi_best_reads);
+ print_in_box(80, 0, 0, " Unique mapping : %s", context->config.report_multi_mapping_reads?"no":"yes");
if(context->config.max_insertion_at_junctions)
- print_in_box(80, 0, 0, "Insertions at junc : %d", context->config.max_insertion_at_junctions);
+ print_in_box(80, 0, 0, " Insertions at junc : %d", context->config.max_insertion_at_junctions);
if(context->config.read_group_id[0])
- print_in_box(80, 0, 0, " Read group name : %s", context->config.read_group_id);
+ print_in_box(80, 0, 0, " Read group name : %s", context->config.read_group_id);
print_in_box(80, 0, 1, "");
print_in_box(80, 2, 1, "http://subread.sourceforge.net/");
@@ -3657,6 +3732,8 @@ int load_global_context(global_context_t * context)
context->config.convert_color_to_base=0;
}
+
+
if(context->input_reads.is_paired_end_reads)
{
if(core_geinput_open(context, &context->input_reads.second_read_file, 2,1))
@@ -3664,7 +3741,6 @@ int load_global_context(global_context_t * context)
//sublog_printf(SUBLOG_STAGE_RELEASED, SUBLOG_LEVEL_ERROR,"Unable to open '%s' as input. Please check if it exists, you have the permission to read it, and it is in the correct format.\n", context->config.second_read_file);
return -1;
}
-
context -> config.max_vote_combinations = 3;
context -> config.multi_best_reads = 3;
context -> config.max_vote_simples = 64;
@@ -3680,24 +3756,14 @@ int load_global_context(global_context_t * context)
context -> config.max_vote_simples = max(context -> config.max_vote_simples , context -> config.reported_multi_best_reads);
context -> config.max_vote_combinations = max(context -> config.max_vote_combinations , context -> config.reported_multi_best_reads);
- if(context->config.reads_per_chunk > 384*1024*1024){
- if(context->input_reads.is_paired_end_reads)
- context->config.reads_per_chunk = 512*1024*1024*min(40,max(0.01,context->config.memory_use_multiplex));
- else
- context->config.reads_per_chunk = 1024*1024*1024*min(40,max(0.01,context->config.memory_use_multiplex));
- }else{
- if(context->input_reads.is_paired_end_reads) context->config.reads_per_chunk /= 2;
- if(context->config.multi_best_reads>1) context->config.reads_per_chunk /= context->config.multi_best_reads;
- //#warning "COMMENT NEXT LINE!!!!!!"
- // context->config.reads_per_chunk /= 16;
- }
+ if(context->input_reads.is_paired_end_reads) context->config.reads_per_chunk /= 2;
+ if(context->config.multi_best_reads>1) context->config.reads_per_chunk /= context->config.multi_best_reads;
struct stat ginp1_stat;
int guess_tested_reads = 0;
stat(context->config.first_read_file , &ginp1_stat);
context->input_reads.first_read_file_size = ginp1_stat.st_size;
-
context -> input_reads.avg_read_length = guess_reads_density_format(context->config.first_read_file , context->config.is_SAM_file_input?1:0, &min_phred_score, &max_phred_score , &guess_tested_reads);
if(context -> input_reads.avg_read_length<0 )context -> input_reads.avg_read_length = 250;
// SUBREADprintf("QR=[%d,%d]; ALEN=%f\n", min_phred_score, max_phred_score, context -> input_reads.avg_read_length);
@@ -4121,8 +4187,8 @@ int chimeric_cigar_parts(global_context_t * global_context, unsigned int sel_pos
{
char * curr_chr, * new_chr;
int curr_offset, new_offset;
- locate_gene_position_max(current_perfect_cursor, &global_context -> chromosome_table, & curr_chr, & curr_offset, 1);
- locate_gene_position_max(jummped_location , &global_context -> chromosome_table, & new_chr, & new_offset, 1);
+ locate_gene_position_max(current_perfect_cursor, &global_context -> chromosome_table, & curr_chr, & curr_offset, NULL, NULL, 1);
+ locate_gene_position_max(jummped_location , &global_context -> chromosome_table, & new_chr, & new_offset, NULL, NULL, 1);
assert(curr_chr);
assert(new_chr);
is_chro_jump = (curr_chr != new_chr);
@@ -4246,7 +4312,7 @@ void quick_sort_run(void * arr, int spot_low,int spot_high, int compare (void *
int pivot,j,i;
if(spot_high-spot_low<1) return;
- pivot = spot_low;
+ pivot = (spot_low + spot_high)/2;
i = spot_low;
j = spot_high;
@@ -4274,7 +4340,24 @@ void quick_sort_run(void * arr, int spot_low,int spot_high, int compare (void *
quick_sort_run(arr, i, spot_high, compare, exchange);
}
+void basic_sort_run(void * arr, int start, int items, int compare (void * arr, int l, int r), void exchange(void * arr, int l, int r)){
+ int i, j;
+ for(i=start; i< start + items - 1; i++)
+ {
+ int min_j = i;
+ for(j=i + 1; j< start + items; j++)
+ {
+ if(compare(arr, min_j, j) > 0)
+ min_j = j;
+ }
+ if(i!=min_j)
+ exchange(arr, i, min_j);
+ }
+}
+void basic_sort(void * arr, int items, int compare (void * arr, int l, int r), void exchange(void * arr, int l, int r)){
+ basic_sort_run(arr, 0, items, compare, exchange);
+}
void merge_sort_run(void * arr, int start, int items, int compare (void * arr, int l, int r), void exchange(void * arr, int l, int r), void merge(void * arr, int start, int items, int items2))
@@ -4288,18 +4371,7 @@ void merge_sort_run(void * arr, int start, int items, int compare (void * arr, i
}
else
{
- int i, j;
- for(i=start; i< start + items - 1; i++)
- {
- int min_j = i;
- for(j=i + 1; j< start + items; j++)
- {
- if(compare(arr, min_j, j) > 0)
- min_j = j;
- }
- if(i!=min_j)
- exchange(arr, i, min_j);
- }
+ basic_sort_run(arr, start, items, compare, exchange);
}
}
void merge_sort(void * arr, int arr_size, int compare (void * arr, int l, int r), void exchange(void * arr, int l, int r), void merge(void * arr, int start, int items, int items2))
diff --git a/src/core.h b/src/core.h
index 28b570b..e521226 100644
--- a/src/core.h
+++ b/src/core.h
@@ -371,7 +371,7 @@ typedef struct{
short best_second_diff_bases;
short realign_flags;
short final_quality;
- short hamming_matched;
+ short chromosomal_length;
} realignment_result_t;
@@ -623,6 +623,7 @@ int chimeric_cigar_parts(global_context_t * global_context , unsigned int sel_po
void warning_file_limit();
void quick_sort(void * arr,int arr_size, int compare (void * arr, int l, int r), void exchange(void * arr, int l, int r));
+void basic_sort(void * arr,int arr_size, int compare (void * arr, int l, int r), void exchange(void * arr, int l, int r));
// L_Minus_R should return -1, 0 or 1 when L<R, L==R or L>R.
// The result is from Small to Large.
diff --git a/src/coverage_calc.c b/src/coverage_calc.c
index f4344a8..eb84d2e 100644
--- a/src/coverage_calc.c
+++ b/src/coverage_calc.c
@@ -185,7 +185,7 @@ int covCalc()
} else {
SUBREADprintf("%s:%s %u [%s] :: %u-%u\n", line_buffer , chro, pos, cigar_str, Staring_Points[x1], Staring_Points[x1]+Section_Lengths[x1]);
SUBREADprintf("Read %s overhangs the boundary of chromosome %s (%u >= %u)\n", line_buffer, chro, x2, chrlen);
- exit(-1);
+ //exit(-1);
}
}
}
@@ -279,7 +279,7 @@ int cov_calc_main(int argc, char ** argv)
return 0;
}
- int is_bam = is_certainly_bam_file(input_file_name, NULL);
+ int is_bam = is_certainly_bam_file(input_file_name, NULL, NULL);
if(1==is_bam) is_BAM_input = 1;
else if(is_bam < 0)
diff --git a/src/gene-algorithms.c b/src/gene-algorithms.c
index 4e0c12e..2e0ca36 100644
--- a/src/gene-algorithms.c
+++ b/src/gene-algorithms.c
@@ -438,7 +438,7 @@ unsigned int get_gene_linear(int chrono, int offset, const unsigned int offsets
}
-int locate_gene_position_max(unsigned int linear, const gene_offset_t* offsets , char ** chro_name, int * pos, int rl)
+int locate_gene_position_max(unsigned int linear, const gene_offset_t* offsets , char ** chro_name, int * pos, int * head_cut_length, int * tail_cut_length, int rl)
{
int n = 0;
@@ -462,13 +462,27 @@ int locate_gene_position_max(unsigned int linear, const gene_offset_t* offsets ,
//SUBREADprintf("max=%u <= lim=%u : ACCEPTED.\n", rl + linear , offsets->read_offsets[n] + 16);
// the end of the read should not excess the end of the chromosome
- if(rl + linear > offsets->read_offsets[n] + 16) return 1;
+ if(tail_cut_length == NULL){
+ if(rl + linear > offsets->read_offsets[n] + 15 - offsets -> padding) return 1;
+ } else {
+ (*tail_cut_length) = linear + rl - ( offsets->read_offsets[n] + 15 - offsets -> padding);
+ if( (*tail_cut_length) >= rl )return 1;
+ }
if (n==0)
*pos = linear;
else
*pos = linear - offsets->read_offsets[n-1];
+ if( (*pos) < offsets -> padding ) {
+ if(head_cut_length == NULL || (*pos) + rl <= offsets -> padding){
+ return 1;
+ }else{
+ (*head_cut_length) = offsets -> padding - (*pos);
+ (*pos) = offsets -> padding;
+ }
+ }
+
(*pos) -= offsets -> padding;
*chro_name = (char *)offsets->read_names+n*MAX_CHROMOSOME_NAME_LEN;
@@ -483,7 +497,7 @@ int locate_gene_position_max(unsigned int linear, const gene_offset_t* offsets ,
int locate_gene_position(unsigned int linear, const gene_offset_t* offsets , char ** chro_name, int * pos)
{
- return locate_gene_position_max(linear, offsets, chro_name, pos, 0);
+ return locate_gene_position_max(linear, offsets, chro_name, pos, NULL, NULL, 0);
}
#define _index_vote(key) (((unsigned int)key)%GENE_VOTE_TABLE_SIZE)
diff --git a/src/gene-algorithms.h b/src/gene-algorithms.h
index a839966..52facb1 100644
--- a/src/gene-algorithms.h
+++ b/src/gene-algorithms.h
@@ -33,7 +33,7 @@ void destroy_offsets(gene_offset_t* offsets);
// Return 0 if the linear position is in a reasonable range or -1 if it is out of range.
// The pointer to the name of the chromosome is put into chro_name, and the position in this chromosome is in pos.
int locate_gene_position(unsigned int linear, const gene_offset_t* offsets , char ** chro_name, int * pos);
-int locate_gene_position_max(unsigned int linear, const gene_offset_t* offsets , char ** chro_name, int * pos, int rl);
+int locate_gene_position_max(unsigned int linear, const gene_offset_t* offsets , char ** chro_name, int * pos, int * head_cut, int * tail_cut, int rl);
unsigned int linear_gene_position(const gene_offset_t* offsets , char *chro_name, unsigned int chro_pos);
diff --git a/src/gene-value-index.c b/src/gene-value-index.c
index d00d896..810ed04 100644
--- a/src/gene-value-index.c
+++ b/src/gene-value-index.c
@@ -59,10 +59,13 @@ int gvindex_get(gene_value_index_t * index, gehash_data_t offset)
unsigned int offset_byte, offset_bit;
gvindex_baseno2offset_m(offset, index , offset_byte, offset_bit);
- if(offset_byte >= index-> values_bytes)return 'N';
+ if(offset_byte >= index-> values_bytes -1)return 'N';
unsigned int one_base_value = (index->values [offset_byte]) >> (offset_bit);
+
+ //SUBREADprintf("RECV_BASE=%d (%d - %d)\n",one_base_value & 3, offset_byte , offset_bit);
+
return int2base(one_base_value & 3);
}
@@ -92,7 +95,7 @@ int gvindex_match(gene_value_index_t * index, gehash_data_t offset, gehash_key_t
}
-void gvindex_set (gene_value_index_t * index, gehash_data_t offset, gehash_key_t base_values, int padding)
+void gvindex_set(gene_value_index_t * index, gehash_data_t offset, gehash_key_t base_values, int padding)
{
unsigned int offset_byte, offset_bit;
gvindex_baseno2offset(offset, index , &offset_byte, &offset_bit);
@@ -112,6 +115,7 @@ void gvindex_set (gene_value_index_t * index, gehash_data_t offset, gehash_key_t
index->values [offset_byte] &= mask;
index->values [offset_byte] |= ((base_values >> (30 - i*2))&0x03) << (offset_bit);
+ //SUBREADprintf("PUT_BASE=%d (%d - %d)\n", (base_values >> (30 - i*2))&0x03, offset_byte, offset_bit);
offset_bit +=2;
if(offset_bit >=8)
{
@@ -133,7 +137,7 @@ void gvindex_dump(gene_value_index_t * index, const char filename [])
unsigned int useful_bytes, useful_bits;
gvindex_baseno2offset (index -> length+ index -> start_point, index,&useful_bytes,&useful_bits);
- fwrite(index->values, 1, useful_bytes, fp);
+ fwrite(index->values, 1, useful_bytes+1, fp);
fclose(fp);
}
@@ -153,8 +157,8 @@ int gvindex_load(gene_value_index_t * index, const char filename [])
unsigned int useful_bytes, useful_bits;
index -> start_base_offset = index -> start_point - index -> start_point%4;
gvindex_baseno2offset (index -> length+ index -> start_point, index ,&useful_bytes,&useful_bits);
- index -> values = malloc(useful_bytes);
- index -> values_bytes = useful_bytes;
+ index -> values = malloc(useful_bytes+1);
+ index -> values_bytes = useful_bytes+1;
if(!index->values)
{
SUBREADputs(MESSAGE_OUT_OF_MEMORY);
@@ -162,7 +166,7 @@ int gvindex_load(gene_value_index_t * index, const char filename [])
}
- read_length =fread(index->values, 1, useful_bytes, fp);
+ read_length =fread(index->values, 1, useful_bytes+1, fp);
assert(read_length>0);
fclose(fp);
diff --git a/src/global-reassembly.c b/src/global-reassembly.c
index 2a429a3..be22803 100644
--- a/src/global-reassembly.c
+++ b/src/global-reassembly.c
@@ -261,6 +261,7 @@ int GRA_check_config (GRA_global_context_t * global_context)
void GRA_init_context(GRA_global_context_t * global_context)
{
+ char mac_rand[13];
FILE * system_random_fp;
memset(global_context, 0 , sizeof(GRA_global_context_t));
@@ -269,7 +270,9 @@ void GRA_init_context(GRA_global_context_t * global_context)
assert(system_random_fp);
fread(global_context -> random_seeds, 8 * sizeof(int), 1, system_random_fp);
fclose(system_random_fp);
- sprintf(global_context -> tmp_file_prefix, "GRAtmp-%u-%d", global_context -> random_seeds[0], global_context -> system_pid);
+
+ mac_or_rand_str(mac_rand);
+ sprintf(global_context -> tmp_file_prefix, "GRAtmp-%06d-%s", global_context -> system_pid, mac_rand);
// best for simulation = 40
global_context -> shortest_read_allowed = 40;
diff --git a/src/index-builder.c b/src/index-builder.c
index 5e148a9..a932bbe 100644
--- a/src/index-builder.c
+++ b/src/index-builder.c
@@ -359,6 +359,8 @@ int build_gene_index(const char index_prefix [], char ** chro_files, int chro_fi
next_char = geinput_next_char(&ginp);
if(next_char < 0)
{
+ gvindex_set(&value_array_index, offset - (IS_COLOR_SPACE?0:0), array_int_key, padding_around_contigs);
+
if (next_char == -1) status = NEXT_READ;
if (next_char == -2) status = NEXT_FILE;
if (next_char == -3) return 0;
diff --git a/src/input-files.c b/src/input-files.c
index 8db4710..9b704e3 100644
--- a/src/input-files.c
+++ b/src/input-files.c
@@ -32,6 +32,7 @@
#include <assert.h>
#include "input-files.h"
#include "sambam-file.h"
+#include "HelperFunctions.h"
#include "hashtable.h"
#include "seek-zlib.h"
#include "gene-algorithms.h"
@@ -1471,7 +1472,7 @@ void write_read_block_file(FILE *temp_fp , unsigned int read_number, char *read_
int get_known_chromosomes(char * in_SAM_file, chromosome_t * known_chromosomes)
{
int i, is_first_read_PE;
- int is_BAM = is_certainly_bam_file(in_SAM_file, &is_first_read_PE);
+ int is_BAM = is_certainly_bam_file(in_SAM_file, &is_first_read_PE, NULL);
SamBam_FILE * fp = SamBam_fopen(in_SAM_file,is_BAM?SAMBAM_FILE_BAM:SAMBAM_FILE_SAM);
while(1)
@@ -2513,6 +2514,8 @@ int SAM_pairer_read_BAM_block(FILE * fp, int max_read_len, char * inbuff) {
int SAM_pairer_read_SAM_MB( FILE * fp, int max_read_len, char * inbuff ){
int ret = 0;
+ if(feof(fp)) return 0;
+
while(1){
if(ret >= max_read_len - MIN_BAM_BLOCK_SIZE || feof(fp))break;
int rlen = fread(inbuff +ret , 1, max_read_len - MIN_BAM_BLOCK_SIZE - ret , fp);
@@ -2521,14 +2524,14 @@ int SAM_pairer_read_SAM_MB( FILE * fp, int max_read_len, char * inbuff ){
int x1;
for(x1 = 0; x1 < min(200, rlen); x1++)
if(*(inbuff+ret+x1)<8 || *(inbuff+ret+x1)> 127){
- // SUBREADprintf("NOT_SAM_ACTUALLY\n");
+ SUBREADprintf("NOT_SAM_ACTUALLY\n");
return -1;
}
ret += rlen;
}
}
if(!feof(fp)){
- int nch, last=-1;
+ int nch;
while(1){
nch = fgetc(fp);
if(nch < 0 || nch == '\n'){
@@ -2536,11 +2539,11 @@ int SAM_pairer_read_SAM_MB( FILE * fp, int max_read_len, char * inbuff ){
}else{
inbuff[ret++]=nch;
}
- last = nch;
}
}
- if(inbuff[ret] != '\n') inbuff[ret++]='\n';
+ if(inbuff[ret-1] != '\n') inbuff[ret++]='\n';
inbuff[ret] = 0;
+
return ret;
}
@@ -2634,7 +2637,7 @@ int SAM_pairer_fetch_BAM_block(SAM_pairer_context_t * pairer , SAM_pairer_thread
if(thread_context -> need_find_start){
int test_read_bin = SAM_pairer_find_start(pairer, thread_context);
- if(test_read_bin<1){
+ if(test_read_bin<1 && thread_context -> input_buff_BIN_used >= 32 ){
pairer -> is_bad_format = 1;
SUBREADprintf("BIN REMAIN=%d, BAM USED=%d, BIN GENERATED=%d, BAM REMAIN=%d, TEST_READ_BIN=%d\n", remained_BIN, used_BAM, have, thread_context -> input_buff_SBAM_used - thread_context -> input_buff_SBAM_ptr, test_read_bin);
}
@@ -2692,7 +2695,9 @@ int SAM_pairer_get_next_read_BIN( SAM_pairer_context_t * pairer , SAM_pairer_thr
unsigned int bam_signature;
BAM_next_u32(bam_signature);
BAM_next_u32(pairer -> BAM_l_text);
- char * header_txt = malloc(pairer->BAM_l_text);
+ char * header_txt = NULL;
+
+ if(pairer->BAM_l_text>0) header_txt = malloc(pairer->BAM_l_text);
//SUBREADprintf("HEAD_TXT=%d\n", pairer -> BAM_l_text);
for(x1 = 0 ; x1 < pairer -> BAM_l_text; x1++){
@@ -2720,18 +2725,20 @@ int SAM_pairer_get_next_read_BIN( SAM_pairer_context_t * pairer , SAM_pairer_thr
memcpy(header_txt + ref_bin_len, &l_ref, 4);
ref_bin_len += 4;
- assert(ref_bin_len < pairer -> BAM_l_text);
- //SUBREADprintf("%d-th ref : %s [len=%u]\n", x1, ref_name, l_ref);
+ if(ref_bin_len >= pairer -> BAM_l_text){
+ SUBREADprintf("%d-th ref : %s [len=%u], bin_len=%d < %d\n", x1, ref_name, l_ref, ref_bin_len, pairer -> BAM_l_text);
+ assert(ref_bin_len < pairer -> BAM_l_text);
+ }
}
//exit(0);
pairer -> output_header(pairer, thread_context -> thread_id, 0, pairer -> BAM_n_ref , header_txt , ref_bin_len );
- free(header_txt);
+ if(header_txt) free(header_txt);
pairer -> BAM_header_parsed = 1;
- if(pairer -> display_progress)
- SUBREADprintf("\nThe header was parsed: %d reference sequences were found.\nScanning the input file.\n", pairer -> BAM_n_ref);
+ //if(pairer -> display_progress)
+ // SUBREADprintf("\nThe header was parsed: %d reference sequences were found.\nScanning the input file.\n", pairer -> BAM_n_ref);
SAM_pairer_fetch_BAM_block(pairer, thread_context);
//SUBREADprintf("HEAD_FINISHED, BAD=%d\n", pairer -> is_bad_format);
@@ -4412,8 +4419,12 @@ int SAM_nosort_decompress_next_block(SAM_pairer_context_t * pairer){
#define NOSORT_SAM_next_line {NOSORT_SAM_eof = fgets(line_ptr, NOSORT_SBAM_BUFF_SIZE, pairer -> input_fp);}
-#define NOSORT_REFILL_LOWBAR 6000*1000
-#define NOSORT_REFILL_HIGHBAR 18000*1000
+#if FEATURECOUNTS_BUFFER_SIZE < ( 12*1024*1024 )
+#error "FEATURECOUNTS_BUFFER_SIZE MUST BE GREATER THAN 12MB!!"
+#endif
+
+#define NOSORT_REFILL_LOWBAR ( 3 * 1024 * 1024 )
+#define NOSORT_REFILL_HIGHBAR ( 6 * 1024 * 1024 )
void SAM_nosort_run_once(SAM_pairer_context_t * pairer){
int x1;
@@ -4630,7 +4641,7 @@ void SAM_nosort_run_once(SAM_pairer_context_t * pairer){
}
record_len = strlen(line_ptr);
- // SUBREADprintf("1CHR=%c, ECHR=%d , RL=%d, RINS=%d, USED=%d\n", line_ptr[0], line_ptr[record_len - 1], record_len, this_thread -> reads_in_SBAM, this_thread -> input_buff_SBAM_used);
+ // SUBREADprintf("1CHR=%c, ECHR=%d , RL=%d, RINS=%d, USED=%d, SIZE=%d\n", line_ptr[0], line_ptr[record_len - 1], record_len, this_thread -> reads_in_SBAM, this_thread -> input_buff_SBAM_used, pairer -> input_buff_SBAM_size);
memcpy(this_thread -> input_buff_SBAM + this_thread -> input_buff_SBAM_used , line_ptr, record_len);
this_thread -> input_buff_SBAM_used += record_len;
this_thread -> reads_in_SBAM ++;
@@ -4678,13 +4689,14 @@ int SAM_pairer_run( SAM_pairer_context_t * pairer){
int sort_SAM_create(SAM_sort_writer * writer, char * output_file, char * tmp_path)
{
- char tmp_fname[MAX_FILE_NAME_LENGTH+40];
+ char tmp_fname[MAX_FILE_NAME_LENGTH+40], mac_rand[13];
memset(writer, 0, sizeof(SAM_sort_writer));
old_sig_TERM = signal (SIGTERM, SAM_SORT_SIGINT_hook);
old_sig_INT = signal (SIGINT, SAM_SORT_SIGINT_hook);
- sprintf(writer -> tmp_path, "%s/temp-sort-%06u-%08X-", tmp_path, getpid(), rand());
+ mac_or_rand_str(mac_rand);
+ sprintf(writer -> tmp_path, "%s/temp-sort-%06u-%s-", tmp_path, getpid(), mac_rand);
_SAMSORT_SNP_delete_temp_prefix = writer -> tmp_path;
sprintf(tmp_fname, "%s%s", writer -> tmp_path, "headers.txt");
@@ -5249,10 +5261,10 @@ int is_SAM_unsorted(char * SAM_line, char * tmp_read_name, short * tmp_flag, uns
return 0;
}
-int is_certainly_bam_file(char * fname, int * is_first_read_PE)
+int is_certainly_bam_file(char * fname, int * is_first_read_PE, long long * SAMBAM_header_size)
{
- int read_type = probe_file_type(fname, is_first_read_PE);
+ int read_type = probe_file_type_EX(fname, is_first_read_PE, SAMBAM_header_size);
if(read_type == FILE_TYPE_NONEXIST || read_type == FILE_TYPE_EMPTY || read_type == FILE_TYPE_UNKNOWN)
return -1;
if(read_type == FILE_TYPE_BAM)
@@ -5458,6 +5470,10 @@ int probe_file_type_fast(char * fname){
}
int probe_file_type(char * fname, int * is_first_read_PE)
{
+ return probe_file_type_EX(fname, is_first_read_PE, NULL);
+}
+int probe_file_type_EX(char * fname, int * is_first_read_PE, long long * SAMBAM_header_length)
+{
FILE * fp = f_subr_open(fname, "rb");
if(!fp) return FILE_TYPE_NONEXIST;
@@ -5590,20 +5606,22 @@ int probe_file_type(char * fname, int * is_first_read_PE)
if(fp)fclose(fp);
- //SUBREADprintf("RET=%d, FIRSTPE=%p\n" , ret, is_first_read_PE);
+ //SUBREADprintf("RET=%d, FIRSTPE=%p, SAMLEN=%p\n" , ret, is_first_read_PE, SAMBAM_header_length);
if(FILE_TYPE_BAM == ret || FILE_TYPE_SAM == ret)
- if(is_first_read_PE)
+ if(is_first_read_PE || SAMBAM_header_length)
{
SamBam_FILE * tpfp = SamBam_fopen(fname, (FILE_TYPE_BAM == ret)?SAMBAM_FILE_BAM:SAMBAM_FILE_SAM);
while(1)
{
char * tbr = SamBam_fgets(tpfp, test_buf, 4999, 0);
- if( tpfp -> is_paired_end >= 10)
+ if( is_first_read_PE && tpfp -> is_paired_end >= 10)
(*is_first_read_PE) = tpfp -> is_paired_end - 10;
if(tbr == NULL)break;
if(tbr[0]=='@') continue;
break;
}
+
+ if( SAMBAM_header_length) (*SAMBAM_header_length) = tpfp -> header_length;
SamBam_fclose(tpfp);
}
diff --git a/src/input-files.h b/src/input-files.h
index 7ddf50d..f98ff3b 100644
--- a/src/input-files.h
+++ b/src/input-files.h
@@ -49,6 +49,8 @@
#define FILE_TYPE_EMPTY 999990
#define FILE_TYPE_NONEXIST 999999
+#define FEATURECOUNTS_BUFFER_SIZE (1024*1024*12)
+
#include <stdlib.h>
@@ -274,7 +276,7 @@ void colorread2base(char * read_buffer, int read_len);
int warning_file_type(char * fname, int expected_type);
char color2char(char clr, char c1);
-int is_certainly_bam_file(char * fname, int * is_firstread_PE);
+int is_certainly_bam_file(char * fname, int * is_firstread_PE, long long * SAMBAM_header_length);
unsigned long long int sort_SAM_hash(char * str);
@@ -286,6 +288,7 @@ int probe_file_type_fast(char * fname);
void geinput_seek(gene_input_t * input, gene_inputfile_position_t * pos);
void geinput_tell(gene_input_t * input, gene_inputfile_position_t * pos);
unsigned long long geinput_file_offset( gene_input_t * input);
+int probe_file_type_EX(char * fname, int * is_first_read_PE, long long * SAMBAM_header_length);
int SAM_pairer_create(SAM_pairer_context_t * pairer, int all_threads, int bin_buff_size_per_thread, int BAM_input, int is_Tiny_Mode, int is_single_end_mode, int force_do_not_sort, int display_progress, char * in_file, void (* reset_output_function) (void * pairer), int (* output_header_function) (void * pairer, int thread_no, int is_text, unsigned int items, char * bin, unsigned int bin_len), int (* output_function) (void * pairer, int thread_no, char * rname, char * bin1, char * bin2), [...]
diff --git a/src/makefile.version b/src/makefile.version
index f05fbf3..177d75d 100644
--- a/src/makefile.version
+++ b/src/makefile.version
@@ -1,3 +1,7 @@
-SUBREAD_VERSION="1.5.0-p1"
+SUBREAD_VERSION_BASE=1.5.0-p2
+SUBREAD_VERSION_DATE=$(SUBREAD_VERSION_BASE)-$(shell date +"%d%b%Y")
+SUBREAD_VERSION="$(SUBREAD_VERSION_DATE)"
+SUBREAD_VERSION="$(SUBREAD_VERSION_BASE)"
+
STATIC_MAKE=
#STATIC_MAKE= -static
diff --git a/src/propmapped.c b/src/propmapped.c
index 0c7ac4e..f5d16c6 100644
--- a/src/propmapped.c
+++ b/src/propmapped.c
@@ -35,7 +35,7 @@
#include "sambam-file.h"
#include "input-files.h"
#include "gene-algorithms.h"
-
+#include "HelperFunctions.h"
typedef struct {
@@ -279,8 +279,10 @@ void add_read_flags(propMapped_context * context, FILE * fp, char * read_name, u
int init_PE_sambam(propMapped_context * context)
{
+ char mac_rand[13];
+ mac_or_rand_str(mac_rand);
srand(time(NULL));
- sprintf(context->temp_file_prefix, "prpm-temp-sum-%06u-%05u", getpid(), rand());
+ sprintf(context->temp_file_prefix, "prpm-temp-sum-%06u-%s", getpid(), mac_rand);
_PROPMAPPED_delete_tmp_prefix = context -> temp_file_prefix;
signal (SIGTERM, PROPMAPPED_SIGINT_hook);
@@ -469,7 +471,7 @@ int propmapped(int argc, char ** argv)
return 0;
}
- int is_bam = is_certainly_bam_file(context -> input_file_name, NULL);
+ int is_bam = is_certainly_bam_file(context -> input_file_name, NULL, NULL);
if(1==is_bam){
context -> is_BAM_input = 1;
diff --git a/src/read-repair.c b/src/read-repair.c
index f39d2dc..8b70d82 100644
--- a/src/read-repair.c
+++ b/src/read-repair.c
@@ -6,6 +6,7 @@
#include "subread.h"
#include "core.h"
#include "input-files.h"
+#include "HelperFunctions.h"
void print_usage_pairer(char * cmd){
SUBREADprintf("\nrepair Version %s\n\n", SUBREAD_VERSION);
@@ -107,9 +108,10 @@ int main_read_repair(int argc, char ** argv)
STANDALONE_exit(-1);
}
- srand( (unsigned int) (miltime()*1000));
+ char mac_rand[13];
+ mac_or_rand_str(mac_rand);
- sprintf(rand_prefix, "fsbm-p%06d-%04X%04X%04X", getpid(), rand()&0xffff, rand()&0xffff,rand()&0xffff);
+ sprintf(rand_prefix, "fsbm-p%06d-%s", getpid(), mac_rand);
SAM_pairer_context_t pairer;
SAM_pairer_writer_main_t writer_main;
diff --git a/src/readSummary.c b/src/readSummary.c
index 3f58a85..a1bbc6d 100644
--- a/src/readSummary.c
+++ b/src/readSummary.c
@@ -122,12 +122,6 @@ typedef struct {
unsigned int chunk_read_ptr;
pthread_t thread_object;
- char * input_buffer;
- unsigned int input_buffer_remainder;
- unsigned int input_buffer_write_ptr;
- pthread_spinlock_t input_buffer_lock;
-
-
unsigned short hits_read_start_base1[MAX_HIT_NUMBER];
unsigned short hits_read_start_base2[MAX_HIT_NUMBER];
@@ -218,7 +212,6 @@ typedef struct {
unsigned short thread_number;
fc_thread_thread_context_t * thread_contexts;
int is_all_finished;
- unsigned int input_buffer_max_size;
int sambam_chro_table_items;
SamBam_Reference_Info * sambam_chro_table;
pthread_spinlock_t sambam_chro_table_lock;
@@ -227,6 +220,7 @@ typedef struct {
char * debug_command;
char * unistr_buffer_space;
+ long long max_BAM_header_size;
unsigned int unistr_buffer_size;
unsigned int unistr_buffer_used;
HashTable * junction_features_table;
@@ -236,9 +230,12 @@ typedef struct {
HashTable * annot_chro_name_alias_table; // name in annotation file -> alias name
char alias_file_name[300];
char input_file_name[300];
+ char * input_file_short_name;
char raw_input_file_name[300];
char output_file_name[300];
+ char output_file_path[300];
unsigned char ** gene_name_array; // gene_internal_number -> gene_name
+ int input_file_unique;
HashTable * exontable_chro_table; // gene_name -> fc_chromosome_index_info structure (contains chro_number, feature_number, block_start, block_end, etc)
int exontable_nchrs;
@@ -556,7 +553,10 @@ void print_FC_configuration(fc_thread_global_context_t * global_context, char *
if(next_fn == NULL || strlen(next_fn)<1) break;
nfiles++;
- int file_probe = is_certainly_bam_file(next_fn, NULL);
+ long long BAM_header_size = -1;
+ int file_probe = is_certainly_bam_file(next_fn, NULL, &BAM_header_size);
+ //SUBREADprintf(" >>> %s : header=%lld , BSIZE=%lld\n", next_fn,BAM_header_size, global_context -> max_BAM_header_size );
+ if(BAM_header_size>0) global_context -> max_BAM_header_size = max( global_context -> max_BAM_header_size , BAM_header_size + 180000);
if(file_probe==-1) nNonExistFiles++;
if(file_probe == 1) nBAMfiles++;
}
@@ -582,7 +582,7 @@ void print_FC_configuration(fc_thread_global_context_t * global_context, char *
{
next_fn = strtok_r(nfiles==0?sam_used:NULL, ";", &tmp_ptr1);
if(next_fn == NULL || strlen(next_fn)<1) break;
- int is_first_read_PE = 0 , file_probe = is_certainly_bam_file(next_fn, &is_first_read_PE);
+ int is_first_read_PE = 0 , file_probe = is_certainly_bam_file(next_fn, &is_first_read_PE, NULL);
char file_chr = 'S';
if(file_probe == -1) file_chr = '?';
@@ -596,7 +596,8 @@ void print_FC_configuration(fc_thread_global_context_t * global_context, char *
(*n_input_files) = nfiles;
print_in_box(80,0,0,"");
print_in_box(80,0,0," Output file : %s", out);
- print_in_box(80,0,0," Annotations : %s (%s)", annot, is_GTF?"GTF":"SAF");
+ print_in_box(80,0,0," Summary : %s.summary", out);
+ print_in_box(80,0,0," Annotation : %s (%s)", annot, is_GTF?"GTF":"SAF");
if(isReadSummaryReport)
print_in_box(80,0,0," Assignment details : <input_file>.featureCounts");
if(global_context -> do_junction_counting)
@@ -615,7 +616,7 @@ void print_FC_configuration(fc_thread_global_context_t * global_context, char *
print_in_box(80,0,0,"");
}
- print_in_box(80,0,0," Strand specific : %s", global_context->is_strand_checked?(global_context->is_strand_checked==1?"yes":"inversed"):"no");
+ print_in_box(80,0,0," Strand specific : %s", global_context->is_strand_checked?(global_context->is_strand_checked==1?"stranded":"reversely stranded"):"no");
char * multi_mapping_allow_mode = "not counted";
if(global_context->is_multi_mapping_allowed == ALLOW_PRIMARY_MAPPING)
multi_mapping_allow_mode = "primary only";
@@ -634,7 +635,7 @@ void print_FC_configuration(fc_thread_global_context_t * global_context, char *
print_in_box(80,0,0," Read reduction to : %d' end" , global_context -> reduce_5_3_ends_to_one == REDUCE_TO_5_PRIME_END ?5:3);
if(global_context -> is_duplicate_ignored)
print_in_box(80,0,0," Duplicated Reads : ignored");
- print_in_box(80,0,0," Read orientations : %c%c", global_context->is_first_read_reversed?'r':'f', global_context->is_second_read_straight?'f':'r' );
+ //print_in_box(80,0,0," Read orientations : %c%c", global_context->is_first_read_reversed?'r':'f', global_context->is_second_read_straight?'f':'r' );
if(global_context->is_paired_end_mode_assign)
{
@@ -651,6 +652,8 @@ void print_FC_configuration(fc_thread_global_context_t * global_context, char *
SUBREADputs("");
print_in_box(80,1,1,"Running");
print_in_box(80,0,0,"");
+ if( global_context -> max_BAM_header_size > 32 * 1024 * 1024 ){
+ }
if(global_context->annot_chro_name_alias_table)
print_in_box(80,0,0,"%ld chromosome name aliases are loaded.", global_context -> annot_chro_name_alias_table ->numOfElements);
@@ -911,15 +914,35 @@ int load_feature_info(fc_thread_global_context_t *global_context, const char * a
unsigned int chro_name_pos = unistr_cpy(global_context, (char *)seq_name, chro_name_len);
global_context -> longest_chro_name = max(chro_name_len, global_context -> longest_chro_name);
+
+ char * start_ptr = strtok_r(NULL,"\t", &token_temp);
+ char * end_ptr = strtok_r(NULL,"\t", &token_temp);
+
+ if(start_ptr == NULL || end_ptr == NULL){
+ SUBREADprintf("\nWarning: the format on the %d-th line is wrong.\n", lineno);
+ }
+ long long int tv1 = atoll(start_ptr);
+ long long int tv2 = atoll(end_ptr);
+
+ if( isdigit(start_ptr[0]) && isdigit(end_ptr[0]) ){
+ if(strlen(start_ptr) > 10 || strlen(end_ptr) > 10 || tv1 > 0x7fffffff || tv2> 0x7fffffff){
+ SUBREADprintf("\nError: Line %d contains a coordinate greater than 2^31!\n", lineno);
+ return -2;
+ }
+ }else{
+ SUBREADprintf("\nError: Line %d contains a format error. The expected annotation format is SAF.\n", lineno);
+ return -2;
+ }
+
ret_features[xk1].chro_name_pos_delta = chro_name_pos - ret_features[xk1].feature_name_pos;
- ret_features[xk1].start = atoi(strtok_r(NULL,"\t", &token_temp));// start
+ ret_features[xk1].start = atoi( start_ptr );// start
if(ret_features[xk1].start<0 || ret_features[xk1].start>0x7fffffff)
{
ret_features[xk1].start = 0;
print_in_box(80,0,0,"WARNING the %d-th line has a negative chro coordinate.", lineno);
}
- ret_features[xk1].end = atoi(strtok_r(NULL,"\t", &token_temp));//end
+ ret_features[xk1].end = atoi( end_ptr );//end
if(ret_features[xk1].end<0 || ret_features[xk1].end>0x7fffffff)
{
ret_features[xk1].end = 0;
@@ -963,8 +986,26 @@ int load_feature_info(fc_thread_global_context_t *global_context, const char * a
char * feature_type = strtok_r(NULL,"\t", &token_temp);// feature_type
if(strcmp(feature_type, global_context -> feature_name_column)==0)
{
- ret_features[xk1].start = atoi(strtok_r(NULL,"\t", &token_temp));// start
- ret_features[xk1].end = atoi(strtok_r(NULL,"\t", &token_temp));//end
+ char * start_ptr = strtok_r(NULL,"\t", &token_temp);
+ char * end_ptr = strtok_r(NULL,"\t", &token_temp);
+
+ if(start_ptr == NULL || end_ptr == NULL){
+ SUBREADprintf("\nWarning: the format on the %d-th line is wrong.\n", lineno);
+ }
+ long long int tv1 = atoll(start_ptr);
+ long long int tv2 = atoll(end_ptr);
+
+ if( isdigit(start_ptr[0]) && isdigit(end_ptr[0]) ){
+ if(strlen(start_ptr) > 10 || strlen(end_ptr) > 10 || tv1 > 0x7fffffff || tv2> 0x7fffffff){
+ SUBREADprintf("\nError: Line %d contains a coordinate greater than 2^31!\n", lineno);
+ return -2;
+ }
+ }else{
+ SUBREADprintf("\nError: Line %d contains a format error. The expected annotation format is GTF/GFF.\n", lineno);
+ return -2;
+ }
+ ret_features[xk1].start = atoi(start_ptr);// start
+ ret_features[xk1].end = atoi(end_ptr);//end
if(ret_features[xk1].start < 1 || ret_features[xk1].end<1 || ret_features[xk1].start > 0x7fffffff || ret_features[xk1].end > 0x7fffffff || ret_features[xk1].start > ret_features[xk1].end)
SUBREADprintf("\nWarning: the feature on the %d-th line has zero coordinate or zero lengths\n\n", lineno);
@@ -2085,8 +2126,28 @@ void process_line_buffer(fc_thread_global_context_t * global_context, fc_thread_
(*cigar_read_len) = Starting_Read_Points[cigar_sections-1] + Section_Read_Lengths[cigar_sections-1];
- if( (is_junction_read && 1 == global_context->is_split_or_exonic_only) || (2 == global_context->is_split_or_exonic_only && !is_junction_read) || !global_context->is_split_or_exonic_only)
- {
+
+ if(global_context->is_split_or_exonic_only == 1 && !is_junction_read) {
+ skipped_for_exonic ++;
+
+ if(skipped_for_exonic == 1 + global_context -> is_paired_end_mode_assign){
+ if(global_context -> SAM_output_fp)
+ fprintf(global_context -> SAM_output_fp,"%s\tUnassigned_%s\t*\t*\n", read_name, (global_context->is_split_or_exonic_only == 2)?"Hasjunction":"Nonjunction");
+
+ thread_context->read_counters.unassigned_junction_condition ++;
+ return;
+ }
+ }
+
+
+ if(global_context->is_split_or_exonic_only == 2 && is_junction_read) {
+ if(global_context -> SAM_output_fp)
+ fprintf(global_context -> SAM_output_fp,"%s\tUnassigned_%s\t*\t*\n", read_name, (global_context->is_split_or_exonic_only == 2)?"Hasjunction":"Nonjunction");
+ thread_context->read_counters.unassigned_junction_condition ++;
+ return;
+ }
+
+ if(1) {
//#warning "=================== COMMENT THESE 2 LINES ================================"
//for(cigar_section_id = 0; cigar_section_id<cigar_sections; cigar_section_id++)
// SUBREADprintf("BCCC: %llu , sec[%d] %s: %u ~ %u ; secs=%d ; flags=%d ; second=%d\n", read_pos, cigar_section_id , ChroNames[cigar_section_id] , Starting_Chro_Points[cigar_section_id], Section_Lengths[cigar_section_id], cigar_sections, alignment_masks, is_second_read);
@@ -2248,23 +2309,9 @@ void process_line_buffer(fc_thread_global_context_t * global_context, fc_thread_
}
}
}
- }else if(global_context->is_split_or_exonic_only == 1 && !is_junction_read) {
- skipped_for_exonic ++;
- if((is_second_read && skipped_for_exonic == 2) || (!global_context -> is_paired_end_mode_assign) || (alignment_masks & 0x8))
- {
- if(global_context -> SAM_output_fp)
- fprintf(global_context -> SAM_output_fp,"%s\tUnassigned_%s\t*\t*\n", read_name, (global_context->is_split_or_exonic_only == 2)?"Hasjunction":"Nonjunction");
-
- thread_context->read_counters.unassigned_junction_condition ++;
- return;
- }
- }else if(global_context->is_split_or_exonic_only == 2 && is_junction_read) {
- if(global_context -> SAM_output_fp)
- fprintf(global_context -> SAM_output_fp,"%s\tUnassigned_%s\t*\t*\n", read_name, (global_context->is_split_or_exonic_only == 2)?"Hasjunction":"Nonjunction");
- thread_context->read_counters.unassigned_junction_condition ++;
- return;
}
+
if(is_second_read) nhits2 = nhits;
else nhits1 = nhits;
} // loop for is_second_read
@@ -2608,6 +2655,21 @@ void vote_and_add_count(fc_thread_global_context_t * global_context, fc_thread_t
final_feture_names[0]=0;
for(xk1 = 0; xk1 < decision_number; xk1++)
{
+
+ int this_decision_score;
+ if(global_context -> use_overlapping_break_tie) {
+ this_decision_score = decision_total_lengths[xk1];
+ } else {
+ if(global_context -> restricted_no_multi_overlap) {
+ this_decision_score = (read1_used[xk1] || read2_used[xk1])?1:0;
+ } else {
+ this_decision_score = read1_used[xk1] + read2_used[xk1];
+ }
+ }
+
+ // This change was made on 31/MAR/2016
+ if( this_decision_score < maximum_decision_score ) continue ;
+
long tmp_voter_id = decision_table_exon_ids[xk1];
thread_context->count_table[tmp_voter_id] += fixed_fractional_count;
@@ -2783,8 +2845,10 @@ HashTable * load_alias_table(char * fname)
void fc_thread_init_global_context(fc_thread_global_context_t * global_context, unsigned int buffer_size, unsigned short threads, int line_length , int is_PE_data, int min_pe_dist, int max_pe_dist, int is_gene_level, int is_overlap_allowed, int is_strand_checked, char * output_fname, int is_sam_out, int is_both_end_required, int is_chimertc_disallowed, int is_PE_distance_checked, char *feature_name_column, char * gene_id_column, int min_map_qual_score, int is_multi_mapping_allowed, int i [...]
{
+ int x1;
+
memset(global_context, 0, sizeof(fc_thread_global_context_t));
- global_context -> input_buffer_max_size = buffer_size;
+ global_context -> max_BAM_header_size = buffer_size;
global_context -> all_reads = 0;
global_context -> redo = 0;
global_context -> SAM_output_fp = NULL;
@@ -2803,8 +2867,8 @@ void fc_thread_init_global_context(fc_thread_global_context_t * global_context,
global_context -> is_multi_mapping_allowed = is_multi_mapping_allowed;
global_context -> is_split_or_exonic_only = is_split_or_exonic_only;
global_context -> is_duplicate_ignored = is_duplicate_ignored;
- global_context -> is_first_read_reversed = (pair_orientations[0]=='r');
- global_context -> is_second_read_straight = (pair_orientations[1]=='f');
+ //global_context -> is_first_read_reversed = (pair_orientations[0]=='r');
+ //global_context -> is_second_read_straight = (pair_orientations[1]=='f');
global_context -> reduce_5_3_ends_to_one = reduce_5_3_ends_to_one;
global_context -> do_not_sort = is_not_sort;
@@ -2827,6 +2891,7 @@ void fc_thread_init_global_context(fc_thread_global_context_t * global_context,
global_context -> calculate_overlapping_lengths = (global_context -> fragment_minimum_overlapping > 1) || global_context -> use_overlapping_break_tie;
global_context -> debug_command = debug_command;
global_context -> max_M = max_M;
+ global_context -> max_BAM_header_size = buffer_size;
global_context -> read_counters.unassigned_ambiguous=0;
global_context -> read_counters.unassigned_nofeatures=0;
@@ -2850,6 +2915,19 @@ void fc_thread_init_global_context(fc_thread_global_context_t * global_context,
strcpy(global_context -> feature_name_column,feature_name_column);
strcpy(global_context -> gene_id_column,gene_id_column);
strcpy(global_context -> output_file_name, output_fname);
+ global_context -> output_file_path[0]=0;
+ for( x1 = strlen(output_fname)-1; x1 >= 0; x1 --){
+ if(output_fname[x1]=='/'){
+ memcpy(global_context -> output_file_path, output_fname, x1);
+ global_context -> output_file_path[x1]=0;
+ break;
+ }
+ }
+ if(0 == global_context -> output_file_path[0]){
+ strcpy(global_context -> output_file_path, ".");
+ }
+
+ //SUBREADprintf("OFPP:%s, OFNN:%s\n", global_context -> output_file_path, global_context -> output_file_name);
global_context -> min_paired_end_distance = min_pe_dist;
global_context -> max_paired_end_distance = max_pe_dist;
@@ -2885,16 +2963,25 @@ int fc_thread_start_threads(fc_thread_global_context_t * global_context, int et_
{
char tmp_fname[350], *modified_fname;
int i=0;
- sprintf(tmp_fname, "%s.featureCounts", global_context -> raw_input_file_name);
- modified_fname = tmp_fname;
- while(modified_fname[0]=='/' || modified_fname[0]=='.' || modified_fname[0]=='\\'){
- modified_fname ++;
- }
- while(modified_fname[i]){
- if(modified_fname[i]=='\\' || modified_fname[i]=='/'||modified_fname[i]==' ')modified_fname[i]='.';
- i++;
+ if( global_context -> input_file_unique ){
+ sprintf(tmp_fname, "%s/%s.featureCounts", global_context -> output_file_path, global_context -> input_file_short_name);
+ global_context -> SAM_output_fp = f_subr_open(tmp_fname, "w");
+ //SUBREADprintf("FCSSF=%s\n", tmp_fname);
+ } else {
+ sprintf(tmp_fname, "%s.featureCounts", global_context -> raw_input_file_name);
+ modified_fname = tmp_fname;
+ while(modified_fname[0]=='/' || modified_fname[0]=='.' || modified_fname[0]=='\\'){
+ modified_fname ++;
+ }
+ while(modified_fname[i]){
+ if(modified_fname[i]=='\\' || modified_fname[i]=='/'||modified_fname[i]==' ')modified_fname[i]='.';
+ i++;
+ }
+ char tmp_fname2[350];
+ sprintf(tmp_fname2, "%s/%s", global_context -> output_file_path, modified_fname);
+ global_context -> SAM_output_fp = f_subr_open(tmp_fname2, "w");
+ //SUBREADprintf("FCSSF=%s\n", tmp_fname2);
}
- global_context -> SAM_output_fp = f_subr_open(modified_fname, "w");
if(!global_context -> SAM_output_fp)
{
SUBREADprintf("Unable to create file '%s'; the read assignment details are not written.\n", tmp_fname);
@@ -2924,10 +3011,6 @@ int fc_thread_start_threads(fc_thread_global_context_t * global_context, int et_
for(xk1=0; xk1<global_context -> thread_number; xk1++)
{
// printf("CHRR_MALLOC\n");
- pthread_spin_init(&global_context->thread_contexts[xk1].input_buffer_lock, PTHREAD_PROCESS_PRIVATE);
- global_context -> thread_contexts[xk1].input_buffer_remainder = 0;
- global_context -> thread_contexts[xk1].input_buffer_write_ptr = 0;
- global_context -> thread_contexts[xk1].input_buffer = malloc(global_context -> input_buffer_max_size);
global_context -> thread_contexts[xk1].thread_id = xk1;
global_context -> thread_contexts[xk1].chunk_read_ptr = 0;
global_context -> thread_contexts[xk1].count_table = calloc(sizeof(read_count_type_t), et_exons);
@@ -2972,9 +3055,12 @@ int fc_thread_start_threads(fc_thread_global_context_t * global_context, int et_
}
char rand_prefix[200];
- sprintf(rand_prefix, "./temp-core-%06u-%08X.sam", getpid(), rand());
+ char MAC_or_random[13];
+ mac_or_rand_str(MAC_or_random);
+ sprintf(rand_prefix, "./temp-core-%06u-%s.sam", getpid(), MAC_or_random);
- SAM_pairer_create(&global_context -> read_pairer, global_context -> thread_number , 64, !global_context-> is_SAM_file, 1, !global_context -> is_paired_end_mode_assign, global_context ->is_paired_end_mode_assign && global_context -> do_not_sort ,0, global_context -> input_file_name, process_pairer_reset, process_pairer_header, process_pairer_output, rand_prefix, global_context);
+ //#warning "REMOVE ' * 2 ' FROM NEXT LINE !!!!!!"
+ SAM_pairer_create(&global_context -> read_pairer, global_context -> thread_number , global_context -> max_BAM_header_size/1024/1024+2, !global_context-> is_SAM_file, 1, !global_context -> is_paired_end_mode_assign, global_context ->is_paired_end_mode_assign && global_context -> do_not_sort ,0, global_context -> input_file_name, process_pairer_reset, process_pairer_header, process_pairer_output, rand_prefix, global_context);
SAM_pairer_set_unsorted_notification(&global_context -> read_pairer, pairer_unsorted_notification);
return 0;
@@ -2995,10 +3081,8 @@ void fc_thread_destroy_thread_context(fc_thread_global_context_t * global_contex
if(global_context -> thread_contexts[xk1].read_coverage_bits)
free(global_context -> thread_contexts[xk1].read_coverage_bits);
free(global_context -> thread_contexts[xk1].count_table);
- free(global_context -> thread_contexts[xk1].input_buffer);
free(global_context -> thread_contexts[xk1].chro_name_buff);
free(global_context -> thread_contexts[xk1].strm_buffer);
- pthread_spin_destroy(&global_context -> thread_contexts[xk1].input_buffer_lock);
if(global_context -> do_junction_counting){
HashTableDestroy(global_context -> thread_contexts[xk1].junction_counting_table);
HashTableDestroy(global_context -> thread_contexts[xk1].splicing_point_table);
@@ -3338,7 +3422,7 @@ static struct option long_options[] =
{"read2pos", required_argument, 0, 0},
{"minOverlap", required_argument, 0, 0},
{"countSplitAlignmentsOnly", no_argument, 0, 0},
- {"countExonicAlignmentsOnly", no_argument, 0, 0},
+ {"countNonSplitAlignmentsOnly", no_argument, 0, 0},
{"debugCommand", required_argument, 0, 0},
{"ignoreDup", no_argument, 0, 0},
{"donotsort", no_argument, 0, 0},
@@ -3358,8 +3442,8 @@ void print_usage()
SUBREADputs("Usage: featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2] ... \n");
SUBREADputs("Required arguments:");
SUBREADputs("");
- SUBREADputs(" -a <string> Name of an annotation file. GTF format by default. See -F ");
- SUBREADputs(" option for more formats.");
+ SUBREADputs(" -a <string> Name of an annotation file. GTF/GFF format by default.");
+ SUBREADputs(" See -F option for more formats.");
SUBREADputs("");
SUBREADputs(" -o <string> Name of the output file including read counts. A separate ");
SUBREADputs(" file including summary statistics of counting results is ");
@@ -3375,9 +3459,9 @@ void print_usage()
SUBREADputs(" with those used in BAM/SAM input, if they are different. ");
SUBREADputs(" See Users Guide for file format.");
SUBREADputs("");
- SUBREADputs(" -F <string> Specify format of provided annotation file. Acceptable ");
- SUBREADputs(" formats include `GTF' and `SAF'. `GTF' by default. See ");
- SUBREADputs(" Users Guide for description of SAF format.");
+ SUBREADputs(" -F <string> Specify format of provided annotation file. Acceptable");
+ SUBREADputs(" formats include `GTF/GFF' and `SAF'. `GTF/GFF' by default.");
+ SUBREADputs(" See Users Guide for description of SAF format.");
SUBREADputs("");
SUBREADputs(" -t <string> Specify feature type in GTF annotation. `exon' by ");
SUBREADputs(" default. Features used for read counting will be ");
@@ -3459,9 +3543,9 @@ void print_usage()
SUBREADputs(" CIGAR string containing `N'). An example of split ");
SUBREADputs(" alignments is exon-spanning reads in RNA-seq data.");
SUBREADputs("");
- SUBREADputs(" -p Count fragments (read pairs) instead of individual reads. ");
- SUBREADputs(" For each read pair, its two reads must be adjacent to ");
- SUBREADputs(" each other in BAM/SAM input.");
+ SUBREADputs(" -p If specified, fragments (or templates) will be counted");
+ SUBREADputs(" instead of reads. This option is only applicable for");
+ SUBREADputs(" paired-end reads.");
SUBREADputs("");
SUBREADputs(" -P Check validity of paired-end distance when counting read ");
SUBREADputs(" pairs. Use -d and -D to set thresholds.");
@@ -3473,9 +3557,6 @@ void print_usage()
SUBREADputs(" -B Count read pairs that have both ends successfully aligned ");
SUBREADputs(" only.");
SUBREADputs("");
- SUBREADputs(" -S <ff:fr:rf> Specify orientation of two reads from the same pair, 'fr' ");
- SUBREADputs(" by by default (forward/reverse).");
- SUBREADputs("");
SUBREADputs(" -C Do not count read pairs that have their two ends mapping ");
SUBREADputs(" to different chromosomes or mapping to same chromosome ");
SUBREADputs(" but on different strands.");
@@ -3831,6 +3912,19 @@ void fc_write_final_junctions(fc_thread_global_context_t * global_context, char
HashTableDestroy(merged_splicing_table);
}
+char * get_short_fname(char * lname){
+ char * ret = lname;
+
+ int x1;
+ for(x1 = strlen(lname)-1; x1>=0; x1--){
+ if(lname [x1] == '/'){
+ ret = lname + x1 + 1;
+ break;
+ }
+ }
+ return ret;
+}
+
void sort_bucket_table(fc_thread_global_context_t * global_context);
int readSummary_single_file(fc_thread_global_context_t * global_context, read_count_type_t * column_numbers, int nexons, int * geneid, char ** chr, long * start, long * stop, unsigned char * sorted_strand, char * anno_chr_2ch, char ** anno_chrs, long * anno_chr_head, long * block_end_index, long * block_min_start , long * block_max_end, fc_read_counters * my_read_counter, HashTable * junc_glob_tab, HashTable * splicing_glob_tab);
@@ -3889,7 +3983,7 @@ int readSummary(int argc,char *argv[]){
long *start, *stop;
int *geneid;
- char *nameFeatureTypeColumn, *nameGeneIDColumn,*debug_command, *pair_orientations;
+ char *nameFeatureTypeColumn, *nameGeneIDColumn,*debug_command, *pair_orientations="fr";
long nexons;
@@ -4028,9 +4122,10 @@ int readSummary(int argc,char *argv[]){
else useOverlappingBreakTie = 0;
- if(argc>35)
+ /*if(argc>35) "-S" is depreciated.
pair_orientations = argv[35];
else pair_orientations = "FR";
+ */
if(argc>36)
doJuncCounting = atoi(argv[36]);
@@ -4051,11 +4146,12 @@ int readSummary(int argc,char *argv[]){
- unsigned int buffer_size = 1024*1024*6;
fc_thread_global_context_t global_context;
- fc_thread_init_global_context(& global_context, buffer_size, thread_number, MAX_LINE_LENGTH, isPE, minPEDistance, maxPEDistance,isGeneLevel, isMultiOverlapAllowed, isStrandChecked, (char *)argv[3] , isReadSummaryReport, isBothEndRequired, isChimericDisallowed, isPEDistChecked, nameFeatureTypeColumn, nameGeneIDColumn, minMappingQualityScore,isMultiMappingAllowed, 0, alias_file_name, cmd_rebuilt, isInputFileResortNeeded, feature_block_size, isCVersion, fiveEndExtension, threeEndExtension [...]
+
+ fc_thread_init_global_context(& global_context, FEATURECOUNTS_BUFFER_SIZE, thread_number, MAX_LINE_LENGTH, isPE, minPEDistance, maxPEDistance,isGeneLevel, isMultiOverlapAllowed, isStrandChecked, (char *)argv[3] , isReadSummaryReport, isBothEndRequired, isChimericDisallowed, isPEDistChecked, nameFeatureTypeColumn, nameGeneIDColumn, minMappingQualityScore,isMultiMappingAllowed, 0, alias_file_name, cmd_rebuilt, isInputFileResortNeeded, feature_block_size, isCVersion, fiveEndExtension, thre [...]
+
if( global_context.is_multi_mapping_allowed != ALLOW_ALL_MULTI_MAPPING && global_context.use_fraction_multi_mapping)
{
@@ -4079,7 +4175,7 @@ int readSummary(int argc,char *argv[]){
print_in_box(84,0,0,"Load annotation file %s %c[0m...", argv[1], CHAR_ESC);
nexons = load_feature_info(&global_context,argv[1], isGTF?FILE_TYPE_GTF:FILE_TYPE_RSUBREAD, &loaded_features);
if(nexons<1){
- SUBREADprintf("Failed to open the annotation file %s, or its format is incorrect, or it contains no '%s' features.\n",argv[1], nameFeatureTypeColumn);
+ if(nexons >= -1) SUBREADprintf("Failed to open the annotation file %s, or its format is incorrect, or it contains no '%s' features.\n",argv[1], nameFeatureTypeColumn);
return -1;
}
@@ -4107,20 +4203,42 @@ int readSummary(int argc,char *argv[]){
}else global_context.fasta_contigs = NULL;
-
global_context.exontable_exons = nexons;
- unsigned int * nreads = (unsigned int *) calloc(nexons,sizeof(int));
-
-
+ unsigned int x1, * nreads = (unsigned int *) calloc(nexons,sizeof(int));
+ char * tmp_pntr = NULL;
+ char * file_list_used = malloc(strlen(argv[2])+1);
+ char * file_list_used2 = malloc(strlen(argv[2])+1);
+ char * is_unique = malloc(strlen(argv[2])+1);
+ strcpy(file_list_used, argv[2]);
+ for(x1 = 0;;x1++){
+ char * test_fn = strtok_r(x1?NULL:file_list_used,";", &tmp_pntr);
+ if(NULL == test_fn) break;
+ char * short_fname = get_short_fname(test_fn);
+ strcpy(file_list_used2, argv[2]);
+
+ is_unique[x1]=1;
+ char * loop_ptr = NULL;
+ int x2;
+ for(x2 = 0;;x2++){
+ char * test_loopfn = strtok_r(x2?NULL:file_list_used2, ";", &loop_ptr);
+ if(NULL == test_loopfn) break;
+ if(x1==x2)continue;
+ char * short_loop_fname = get_short_fname(test_loopfn);
+ if(strcmp(short_loop_fname, short_fname)==0) {
+ is_unique[x1] = 0;
+ break;
+ }
+ }
+ }
+ free(file_list_used2);
- char * tmp_pntr = NULL;
- char * file_list_used = malloc(strlen(argv[2])+1);
+ tmp_pntr = NULL;
strcpy(file_list_used, argv[2]);
char * next_fn = strtok_r(file_list_used,";", &tmp_pntr);
read_count_type_t ** table_columns = calloc( n_input_files , sizeof(read_count_type_t *)), i_files=0;
@@ -4133,7 +4251,7 @@ int readSummary(int argc,char *argv[]){
splicing_global_table_list = calloc(n_input_files, sizeof(HashTable *));
}
- for(;;){
+ for(x1 = 0;;x1++){
int orininal_isPE = global_context.is_paired_end_mode_assign;
if(next_fn==NULL || strlen(next_fn)<1) break;
@@ -4143,6 +4261,9 @@ int readSummary(int argc,char *argv[]){
strcpy(global_context.input_file_name, next_fn);
strcpy(global_context.raw_input_file_name, next_fn);
+ global_context.input_file_unique = is_unique[x1];
+ global_context.input_file_short_name = get_short_fname(next_fn);
+ //SUBREADprintf("UNQQ=%d ; FNAME=%s\n", global_context.input_file_unique, global_context.input_file_name);
global_context.redo=0;
@@ -4338,7 +4459,7 @@ int readSummary_single_file(fc_thread_global_context_t * global_context, read_co
if(strcmp( global_context->input_file_name,"STDIN")!=0)
{
- int file_probe = is_certainly_bam_file(global_context->input_file_name, &is_first_read_PE);
+ int file_probe = is_certainly_bam_file(global_context->input_file_name, &is_first_read_PE, NULL);
global_context -> is_paired_end_input_file = is_first_read_PE;
// a Singel-end SAM/BAM file cannot be assigned as a PE SAM/BAM file;
@@ -4522,6 +4643,7 @@ int feature_count_main(int argc, char ** argv)
switch(c)
{
case 'S':
+ /*
if(strlen(optarg)!=2 || (strcmp(optarg, "ff")!=0 && strcmp(optarg, "rf")!=0 && strcmp(optarg, "fr")!=0)){
SUBREADprintf("The order parameter can only be ff, fr or rf.\n");
print_usage();
@@ -4530,6 +4652,8 @@ int feature_count_main(int argc, char ** argv)
Pair_Orientations[0]=(optarg[0]=='r'?'r':'f');
Pair_Orientations[1]=(optarg[1]=='f'?'f':'r');
Pair_Orientations[2]=0;
+ */
+ SUBREADprintf("The \"-S\" option has been depreciated.\n");
break;
case 'G':
@@ -4707,7 +4831,7 @@ int feature_count_main(int argc, char ** argv)
{
is_Restrictedly_No_Overlap = 1;
}
- if(strcmp("countExonicAlignmentsOnly", long_options[option_index].name)==0)
+ if(strcmp("countNonSplitAlignmentsOnly", long_options[option_index].name)==0)
{
is_Split_or_Exonic_Only = 2;
}
diff --git a/src/removeDupReads.c b/src/removeDupReads.c
index c03e77b..da5e10c 100644
--- a/src/removeDupReads.c
+++ b/src/removeDupReads.c
@@ -25,6 +25,7 @@
#include <stdlib.h>
#include <string.h>
#include "subread.h"
+#include "HelperFunctions.h"
#include "removeDupReads.h"
#include "input-files.h"
extern unsigned int BASE_BLOCK_LENGTH;
@@ -199,7 +200,10 @@ int repeated_read_removal(char * in_SAM_file, int threshold, char * out_SAM_file
memcpy(rand48_seed, &start_time, 6);
seed48(rand48_seed);
- sprintf(temp_file_prefix, "%s/temp-delrep-%06u-%08lX-", temp_location==NULL?".":temp_location, getpid(), lrand48());
+
+ char mac_rand[13];
+ mac_or_rand_str(mac_rand);
+ sprintf(temp_file_prefix, "%s/temp-delrep-%06u-%s-", temp_location==NULL?".":temp_location, getpid(), mac_rand);
if(break_SAM_file(in_SAM_file, 0, temp_file_prefix, &real_read_count, NULL, known_chromosomes, 0 /* This 0 means that the sequence/quality/cigar fields are not needed in the temp files*/, 0, NULL, NULL, NULL, NULL, NULL)) return -1;
diff --git a/src/sambam-file.c b/src/sambam-file.c
index b8339d4..c38250c 100644
--- a/src/sambam-file.c
+++ b/src/sambam-file.c
@@ -149,6 +149,7 @@ SamBam_FILE * SamBam_fopen(char * fname , int file_type)
SB_RINC(ret, 4);
ret -> bam_file_next_section_start = ret -> input_binary_stream_read_ptr + l_text;
+ ret -> header_length = ret -> bam_file_next_section_start;
}
return ret;
}
@@ -268,6 +269,10 @@ char * SamBam_fgets(SamBam_FILE * fp, char * buff , int buff_len, int seq_needed
}
}
}
+
+ if(buff[0] == '@'){
+ fp -> header_length = ftello(fp->os_file) + strlenbuff+1;
+ }
return ret;
}
}
@@ -312,6 +317,7 @@ char * SamBam_fgets(SamBam_FILE * fp, char * buff , int buff_len, int seq_needed
{
SamBam_read_ref_info(fp);
fp -> bam_file_stage = BAM_FILE_STAGE_ALIGNMENT;
+ fp -> header_length = fp-> input_binary_stream_read_ptr;
}
return buff;
}
diff --git a/src/sambam-file.h b/src/sambam-file.h
index e53a183..d8193cb 100644
--- a/src/sambam-file.h
+++ b/src/sambam-file.h
@@ -88,6 +88,7 @@ typedef struct
unsigned long long input_binary_stream_read_ptr;
unsigned long long input_binary_stream_write_ptr;
unsigned long long input_binary_stream_buffer_start_ptr;
+ unsigned long long header_length;
SamBam_Reference_Info * bam_chro_table;
int bam_chro_table_size;
diff --git a/src/subread.h b/src/subread.h
index 368a7bf..8b3d959 100644
--- a/src/subread.h
+++ b/src/subread.h
@@ -211,7 +211,7 @@ typedef short gene_vote_number_t;
//#define base2int(c) ("\x0\x0\x2\x0\x0\x0\x1\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3\x3"[(c)-'A'])
//#define int2base(c) ((c)==1?'G':((c)==0?'A':((c)==2?'C':'T')))
//#define int2base(c) ("AGCT"[(c)])
-#define int2base(c) (1413695297 >> (8*(c))&0xff);
+#define int2base(c) (1413695297 >> (8*(c))&0xff)
#define color2int(c) ((c) - '0')
#define int2color(c) ("0123"[(c)])
#define remove_backslash(str) { int xxxa=0; while(str[xxxa]){ if(str[xxxa]=='/'){str[xxxa]='\0'; break;} xxxa++;} }
diff --git a/src/subtools.c b/src/subtools.c
index 6b19bbb..17fa0ca 100644
--- a/src/subtools.c
+++ b/src/subtools.c
@@ -127,10 +127,11 @@ int main(int argc, char ** argv)
if(sort_needed)
{
- char fline[3000], temp_file_name[300];
+ char fline[3000], temp_file_name[300], mac_rand[13];
SAM_sort_writer writer;
- sprintf(temp_file_name, "./temp-subt-%06u-%08X.sam", getpid(), rand());
+ mac_or_rand_str(mac_rand);
+ sprintf(temp_file_name, "./temp-subt-%06u-%s.sam", getpid(), mac_rand);
int ret =sort_SAM_create(&writer, temp_file_name, ".");
if(ret)
{
diff --git a/test/featureCounts/data/corner-EXON-ONLY.ora b/test/featureCounts/data/corner-EXON-ONLY.ora
new file mode 100644
index 0000000..ac3835f
--- /dev/null
+++ b/test/featureCounts/data/corner-EXON-ONLY.ora
@@ -0,0 +1,9 @@
+# Program
+Geneid Chr Start End corner-JUNC.sam
+simu_gene1 0
+simu_gene2 1
+simu_gene3 0
+simu_gene4 1
+simu_gene5 0
+simu_gene6 0
+simu_gene7 0
diff --git a/test/featureCounts/data/test-minimum-dup.ora.summary b/test/featureCounts/data/test-minimum-dup.ora.summary
deleted file mode 100644
index 6f8f328..0000000
--- a/test/featureCounts/data/test-minimum-dup.ora.summary
+++ /dev/null
@@ -1,12 +0,0 @@
-Status test-chrname.sam
-Assigned 139
-Unassigned_Ambiguity 1
-Unassigned_MultiMapping 0
-Unassigned_NoFeatures 158
-Unassigned_Unmapped 0
-Unassigned_MappingQuality 0
-Unassigned_FragementLength 0
-Unassigned_Chimera 0
-Unassigned_Secondary 0
-Unassigned_Nonjunction 0
-Unassigned_Duplicate 302
diff --git a/test/featureCounts/test_corner_cases.sh b/test/featureCounts/test_corner_cases.sh
index 0a7fc59..12c9843 100644
--- a/test/featureCounts/test_corner_cases.sh
+++ b/test/featureCounts/test_corner_cases.sh
@@ -58,5 +58,6 @@ $SH_CMD data/compare.sh data/test-minimum.sam data/test-minimum-UNSTR.ora data/t
# test 5' and 3' end extension
$SH_CMD data/compare.sh data/test-chrname.sam data/test-minimum-dup.ora data/test-minimum.GTF " -p --ignoreDup " "Ignoring duplicate fragments"
$SH_CMD data/compare.sh data/corner-JUNC.sam data/corner-JUNC-ONLY.ora data/test-minimum.GTF "--countSplitAlignmentsOnly -O -f " "Junction reads only" FL
+$SH_CMD data/compare.sh data/corner-JUNC.sam data/corner-EXON-ONLY.ora data/test-minimum.GTF "--countNonSplitAlignmentsOnly -p " "Exonic reads only"
echo
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/subread.git
More information about the debian-med-commit
mailing list