Commit a01f6ca1 authored by Charles Plessy's avatar Charles Plessy

Imported Upstream version 0.5.10

parent dbe76989
......@@ -163,7 +163,7 @@ void bns_destroy(bntseq_t *bns)
}
}
void bns_fasta2bntseq(gzFile fp_fa, const char *prefix)
int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix)
{
kseq_t *seq;
char name[1024];
......@@ -172,6 +172,7 @@ void bns_fasta2bntseq(gzFile fp_fa, const char *prefix)
int l_buf;
unsigned char buf[0x10000];
int32_t m_seqs, m_holes, l, i;
int64_t ret = -1;
FILE *fp;
// initialization
......@@ -235,6 +236,7 @@ void bns_fasta2bntseq(gzFile fp_fa, const char *prefix)
bns->l_pac += seq->seq.l;
}
xassert(bns->l_pac, "zero length sequence.");
ret = bns->l_pac;
{ // finalize .pac file
ubyte_t ct;
fwrite(buf, 1, (l_buf>>2) + ((l_buf&3) == 0? 0 : 1), fp);
......@@ -251,6 +253,7 @@ void bns_fasta2bntseq(gzFile fp_fa, const char *prefix)
bns_dump(bns, prefix);
bns_destroy(bns);
kseq_destroy(seq);
return ret;
}
int bwa_fa2pac(int argc, char *argv[])
......
......@@ -70,7 +70,7 @@ extern "C" {
bntseq_t *bns_restore(const char *prefix);
bntseq_t *bns_restore_core(const char *ann_filename, const char* amb_filename, const char* pac_filename);
void bns_destroy(bntseq_t *bns);
void bns_fasta2bntseq(gzFile fp_fa, const char *prefix);
int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix);
int bns_coor_pac2real(const bntseq_t *bns, int64_t pac_coor, int len, int32_t *real_seq);
#ifdef __cplusplus
......
This diff is collapsed.
......@@ -341,7 +341,8 @@ int bwa_cal_pac_pos_pe(const char *prefix, bwt_t *const _bwt[2], int n_seqs, bwa
&& (p[1]->type == BWA_TYPE_UNIQUE || p[1]->type == BWA_TYPE_REPEAT))
{ // only when both ends mapped
uint64_t x;
int j, k, n_occ[2];
int j, k;
long long n_occ[2];
for (j = 0; j < 2; ++j) {
n_occ[j] = 0;
for (k = 0; k < d->aln[j].n; ++k)
......
......@@ -437,15 +437,15 @@ void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, in
if (mate->strand) flag |= SAM_FMR;
} else flag |= SAM_FMU;
}
printf("%s\t%d\t%s\t", p->name, flag, bns->anns[seqid].name);
printf("%d\t%d\t", (int)(p->pos - bns->anns[seqid].offset + 1), p->mapQ);
err_printf("%s\t%d\t%s\t", p->name, flag, bns->anns[seqid].name);
err_printf("%d\t%d\t", (int)(p->pos - bns->anns[seqid].offset + 1), p->mapQ);
// print CIGAR
if (p->cigar) {
for (j = 0; j != p->n_cigar; ++j)
printf("%d%c", __cigar_len(p->cigar[j]), "MIDS"[__cigar_op(p->cigar[j])]);
} else if (p->type == BWA_TYPE_NO_MATCH) printf("*");
else printf("%dM", p->len);
err_printf("%d%c", __cigar_len(p->cigar[j]), "MIDS"[__cigar_op(p->cigar[j])]);
} else if (p->type == BWA_TYPE_NO_MATCH) err_printf("*");
else err_printf("%dM", p->len);
// print mate coordinate
if (mate && mate->type != BWA_TYPE_NO_MATCH) {
......@@ -454,12 +454,12 @@ void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, in
am = mate->seQ < p->seQ? mate->seQ : p->seQ; // smaller single-end mapping quality
// redundant calculation here, but should not matter too much
m_is_N = bns_coor_pac2real(bns, mate->pos, mate->len, &m_seqid);
printf("\t%s\t", (seqid == m_seqid)? "=" : bns->anns[m_seqid].name);
err_printf("\t%s\t", (seqid == m_seqid)? "=" : bns->anns[m_seqid].name);
isize = (seqid == m_seqid)? pos_5(mate) - pos_5(p) : 0;
if (p->type == BWA_TYPE_NO_MATCH) isize = 0;
printf("%d\t%lld\t", (int)(mate->pos - bns->anns[m_seqid].offset + 1), isize);
} else if (mate) printf("\t=\t%d\t0\t", (int)(p->pos - bns->anns[seqid].offset + 1));
else printf("\t*\t0\t0\t");
err_printf("%d\t%lld\t", (int)(mate->pos - bns->anns[m_seqid].offset + 1), isize);
} else if (mate) err_printf("\t=\t%d\t0\t", (int)(p->pos - bns->anns[seqid].offset + 1));
else err_printf("\t*\t0\t0\t");
// print sequence and quality
if (p->strand == 0)
......@@ -468,42 +468,42 @@ void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, in
putchar('\t');
if (p->qual) {
if (p->strand) seq_reverse(p->len, p->qual, 0); // reverse quality
printf("%s", p->qual);
} else printf("*");
err_printf("%s", p->qual);
} else err_printf("*");
if (bwa_rg_id) printf("\tRG:Z:%s", bwa_rg_id);
if (p->bc[0]) printf("\tBC:Z:%s", p->bc);
if (p->clip_len < p->full_len) printf("\tXC:i:%d", p->clip_len);
if (bwa_rg_id) err_printf("\tRG:Z:%s", bwa_rg_id);
if (p->bc[0]) err_printf("\tBC:Z:%s", p->bc);
if (p->clip_len < p->full_len) err_printf("\tXC:i:%d", p->clip_len);
if (p->type != BWA_TYPE_NO_MATCH) {
int i;
// calculate XT tag
XT = "NURM"[p->type];
if (nn > 10) XT = 'N';
// print tags
printf("\tXT:A:%c\t%s:i:%d", XT, (mode & BWA_MODE_COMPREAD)? "NM" : "CM", p->nm);
if (nn) printf("\tXN:i:%d", nn);
if (mate) printf("\tSM:i:%d\tAM:i:%d", p->seQ, am);
err_printf("\tXT:A:%c\t%s:i:%d", XT, (mode & BWA_MODE_COMPREAD)? "NM" : "CM", p->nm);
if (nn) err_printf("\tXN:i:%d", nn);
if (mate) err_printf("\tSM:i:%d\tAM:i:%d", p->seQ, am);
if (p->type != BWA_TYPE_MATESW) { // X0 and X1 are not available for this type of alignment
printf("\tX0:i:%d", p->c1);
if (p->c1 <= max_top2) printf("\tX1:i:%d", p->c2);
err_printf("\tX0:i:%d", p->c1);
if (p->c1 <= max_top2) err_printf("\tX1:i:%d", p->c2);
}
printf("\tXM:i:%d\tXO:i:%d\tXG:i:%d", p->n_mm, p->n_gapo, p->n_gapo+p->n_gape);
if (p->md) printf("\tMD:Z:%s", p->md);
err_printf("\tXM:i:%d\tXO:i:%d\tXG:i:%d", p->n_mm, p->n_gapo, p->n_gapo+p->n_gape);
if (p->md) err_printf("\tMD:Z:%s", p->md);
// print multiple hits
if (p->n_multi) {
printf("\tXA:Z:");
err_printf("\tXA:Z:");
for (i = 0; i < p->n_multi; ++i) {
bwt_multi1_t *q = p->multi + i;
int k;
j = pos_end_multi(q, p->len) - q->pos;
nn = bns_coor_pac2real(bns, q->pos, j, &seqid);
printf("%s,%c%d,", bns->anns[seqid].name, q->strand? '-' : '+',
err_printf("%s,%c%d,", bns->anns[seqid].name, q->strand? '-' : '+',
(int)(q->pos - bns->anns[seqid].offset + 1));
if (q->cigar) {
for (k = 0; k < q->n_cigar; ++k)
printf("%d%c", __cigar_len(q->cigar[k]), "MIDS"[__cigar_op(q->cigar[k])]);
} else printf("%dM", p->len);
printf(",%d;", q->gap + q->mm);
err_printf("%d%c", __cigar_len(q->cigar[k]), "MIDS"[__cigar_op(q->cigar[k])]);
} else err_printf("%dM", p->len);
err_printf(",%d;", q->gap + q->mm);
}
}
}
......@@ -512,16 +512,16 @@ void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, in
ubyte_t *s = p->strand? p->rseq : p->seq;
int flag = p->extra_flag | SAM_FSU;
if (mate && mate->type == BWA_TYPE_NO_MATCH) flag |= SAM_FMU;
printf("%s\t%d\t*\t0\t0\t*\t*\t0\t0\t", p->name, flag);
err_printf("%s\t%d\t*\t0\t0\t*\t*\t0\t0\t", p->name, flag);
for (j = 0; j != p->len; ++j) putchar("ACGTN"[(int)s[j]]);
putchar('\t');
if (p->qual) {
if (p->strand) seq_reverse(p->len, p->qual, 0); // reverse quality
printf("%s", p->qual);
} else printf("*");
if (bwa_rg_id) printf("\tRG:Z:%s", bwa_rg_id);
if (p->bc[0]) printf("\tBC:Z:%s", p->bc);
if (p->clip_len < p->full_len) printf("\tXC:i:%d", p->clip_len);
err_printf("%s", p->qual);
} else err_printf("*");
if (bwa_rg_id) err_printf("\tRG:Z:%s", bwa_rg_id);
if (p->bc[0]) err_printf("\tBC:Z:%s", p->bc);
if (p->clip_len < p->full_len) err_printf("\tXC:i:%d", p->clip_len);
putchar('\n');
}
}
......@@ -541,8 +541,8 @@ void bwa_print_sam_SQ(const bntseq_t *bns)
{
int i;
for (i = 0; i < bns->n_seqs; ++i)
printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[i].name, bns->anns[i].len);
if (bwa_rg_line) printf("%s\n", bwa_rg_line);
err_printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[i].name, bns->anns[i].len);
if (bwa_rg_line) err_printf("%s\n", bwa_rg_line);
}
void bwase_initialize()
......
......@@ -149,14 +149,21 @@ bwa_seq_t *bwa_read_seq(bwa_seqio_t *bs, int n_needed, int *n, int mode, int tri
int n_seqs, l, i, is_comp = mode&BWA_MODE_COMPREAD, is_64 = mode&BWA_MODE_IL13, l_bc = mode>>24;
long n_trimmed = 0, n_tot = 0;
if (l_bc > 15) {
fprintf(stderr, "[%s] the maximum barcode length is 15.\n", __func__);
if (l_bc > BWA_MAX_BCLEN) {
fprintf(stderr, "[%s] the maximum barcode length is %d.\n", __func__, BWA_MAX_BCLEN);
return 0;
}
if (bs->is_bam) return bwa_read_bam(bs, n_needed, n, is_comp, trim_qual); // l_bc has no effect for BAM input
n_seqs = 0;
seqs = (bwa_seq_t*)calloc(n_needed, sizeof(bwa_seq_t));
while ((l = kseq_read(seq)) >= 0) {
if ((mode & BWA_MODE_CFY) && (seq->comment.l != 0)) {
// skip reads that are marked to be filtered by Casava
char *s = index(seq->comment.s, ':');
if (s && *(++s) == 'Y') {
continue;
}
}
if (is_64 && seq->qual.l)
for (i = 0; i < seq->qual.l; ++i) seq->qual.s[i] -= 31;
if (seq->seq.l <= l_bc) continue; // sequence length equals or smaller than the barcode length
......
......@@ -15,7 +15,7 @@ SUBDIRS=
lib:libbwtgen.a
libbwtgen.a:$(OBJS)
$(AR) -cru $@ $(OBJS)
$(AR) -scru $@ $(OBJS)
cleanlocal:
rm -f gmon.out *.o a.out $(PROG) *~ *.a
......
......@@ -13,9 +13,7 @@
#include "utils.h"
#ifdef HAVE_PTHREAD
#define THREAD_BLOCK_SIZE 1024
#include <pthread.h>
static pthread_mutex_t g_seq_lock = PTHREAD_MUTEX_INITIALIZER;
#endif
gap_opt_t *gap_init_opt()
......@@ -98,18 +96,7 @@ void bwa_cal_sa_reg_gap(int tid, bwt_t *const bwt[2], int n_seqs, bwa_seq_t *seq
for (i = 0; i != n_seqs; ++i) {
bwa_seq_t *p = seqs + i;
#ifdef HAVE_PTHREAD
if (opt->n_threads > 1) {
pthread_mutex_lock(&g_seq_lock);
if (p->tid < 0) { // unassigned
int j;
for (j = i; j < n_seqs && j < i + THREAD_BLOCK_SIZE; ++j)
seqs[j].tid = tid;
} else if (p->tid != tid) {
pthread_mutex_unlock(&g_seq_lock);
continue;
}
pthread_mutex_unlock(&g_seq_lock);
}
if (i % opt->n_threads != tid) continue;
#endif
p->sa = 0; p->type = BWA_TYPE_NO_MATCH; p->c1 = p->c2 = 0; p->n_aln = 0; p->aln = 0;
seq[0] = p->seq; seq[1] = p->rseq;
......@@ -189,7 +176,7 @@ void bwa_aln_core(const char *prefix, const char *fn_fa, const gap_opt_t *opt)
}
// core loop
fwrite(opt, sizeof(gap_opt_t), 1, stdout);
err_fwrite(opt, sizeof(gap_opt_t), 1, stdout);
while ((seqs = bwa_read_seq(ks, 0x40000, &n_seqs, opt->mode, opt->trim_qual)) != 0) {
tot_seqs += n_seqs;
t = clock();
......@@ -226,8 +213,8 @@ void bwa_aln_core(const char *prefix, const char *fn_fa, const gap_opt_t *opt)
fprintf(stderr, "[bwa_aln_core] write to the disk... ");
for (i = 0; i < n_seqs; ++i) {
bwa_seq_t *p = seqs + i;
fwrite(&p->n_aln, 4, 1, stdout);
if (p->n_aln) fwrite(p->aln, sizeof(bwt_aln1_t), p->n_aln, stdout);
err_fwrite(&p->n_aln, 4, 1, stdout);
if (p->n_aln) err_fwrite(p->aln, sizeof(bwt_aln1_t), p->n_aln, stdout);
}
fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC); t = clock();
......@@ -246,7 +233,7 @@ int bwa_aln(int argc, char *argv[])
gap_opt_t *opt;
opt = gap_init_opt();
while ((c = getopt(argc, argv, "n:o:e:i:d:l:k:cLR:m:t:NM:O:E:q:f:b012IB:")) >= 0) {
while ((c = getopt(argc, argv, "n:o:e:i:d:l:k:cLR:m:t:NM:O:E:q:f:b012IYB:")) >= 0) {
switch (c) {
case 'n':
if (strstr(optarg, ".")) opt->fnr = atof(optarg), opt->max_diff = -1;
......@@ -274,6 +261,7 @@ int bwa_aln(int argc, char *argv[])
case '1': opt->mode |= BWA_MODE_BAM_READ1; break;
case '2': opt->mode |= BWA_MODE_BAM_READ2; break;
case 'I': opt->mode |= BWA_MODE_IL13; break;
case 'Y': opt->mode |= BWA_MODE_CFY; break;
case 'B': opt->mode |= atoi(optarg) << 24; break;
default: return 1;
}
......@@ -311,6 +299,7 @@ int bwa_aln(int argc, char *argv[])
fprintf(stderr, " -0 use single-end reads only (effective with -b)\n");
fprintf(stderr, " -1 use the 1st read in a pair (effective with -b)\n");
fprintf(stderr, " -2 use the 2nd read in a pair (effective with -b)\n");
fprintf(stderr, " -Y filter Casava-filtered sequences\n");
fprintf(stderr, "\n");
return 1;
}
......
......@@ -22,6 +22,8 @@
#define BWA_AVG_ERR 0.02
#define BWA_MIN_RDLEN 35 // for read trimming
#define BWA_MAX_BCLEN 63 // maximum barcode length; 127 is the maximum
#ifndef bns_pac
#define bns_pac(pac, k) ((pac)[(k)>>2] >> ((~(k)&3)<<1) & 3)
#endif
......@@ -75,7 +77,7 @@ typedef struct {
// for multi-threading only
int tid;
// barcode
char bc[16]; // null terminated; up to 15 bases
char bc[BWA_MAX_BCLEN+1]; // null terminated; up to BWA_MAX_BCLEN bases
// NM and MD tags
uint32_t full_len:20, nm:12;
char *md;
......@@ -84,6 +86,7 @@ typedef struct {
#define BWA_MODE_GAPE 0x01
#define BWA_MODE_COMPREAD 0x02
#define BWA_MODE_LOGGAP 0x04
#define BWA_MODE_CFY 0x08
#define BWA_MODE_NONSTOP 0x10
#define BWA_MODE_BAM 0x20
#define BWA_MODE_BAM_SE 0x40
......
......@@ -57,7 +57,7 @@ static inline void gap_push(gap_stack_t *stack, int a, int i, bwtint_t k, bwtint
p = q->stack + q->n_entries;
p->info = (u_int32_t)score<<21 | a<<20 | i; p->k = k; p->l = l;
p->n_mm = n_mm; p->n_gapo = n_gapo; p->n_gape = n_gape; p->state = state;
if (is_diff) p->last_diff_pos = i;
p->last_diff_pos = is_diff? i : 0;
++(q->n_entries);
++(stack->n_entries);
if (stack->best > score) stack->best = score;
......
......@@ -42,12 +42,13 @@ void bwa_pac_rev_core(const char *fn, const char *fn_rev);
int bwa_index(int argc, char *argv[])
{
char *prefix = 0, *str, *str2, *str3;
int c, algo_type = 3, is_color = 0;
int c, algo_type = 0, is_color = 0;
clock_t t;
int64_t l_pac;
while ((c = getopt(argc, argv, "ca:p:")) >= 0) {
switch (c) {
case 'a':
case 'a': // if -a is not set, algo_type will be determined later
if (strcmp(optarg, "div") == 0) algo_type = 1;
else if (strcmp(optarg, "bwtsw") == 0) algo_type = 2;
else if (strcmp(optarg, "is") == 0) algo_type = 3;
......@@ -79,7 +80,7 @@ int bwa_index(int argc, char *argv[])
gzFile fp = xzopen(argv[optind], "r");
t = clock();
fprintf(stderr, "[bwa_index] Pack FASTA... ");
bns_fasta2bntseq(fp, prefix);
l_pac = bns_fasta2bntseq(fp, prefix);
fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
gzclose(fp);
} else { // color indexing
......@@ -87,7 +88,7 @@ int bwa_index(int argc, char *argv[])
strcat(strcpy(str, prefix), ".nt");
t = clock();
fprintf(stderr, "[bwa_index] Pack nucleotide FASTA... ");
bns_fasta2bntseq(fp, str);
l_pac = bns_fasta2bntseq(fp, str);
fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
gzclose(fp);
{
......@@ -99,6 +100,11 @@ int bwa_index(int argc, char *argv[])
fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC);
}
}
if (l_pac > 0xffffffffu) {
fprintf(stderr, "[%s] BWA only works with reference sequences shorter than 4GB in total. Abort!\n", __func__);
return 1;
}
if (algo_type == 0) algo_type = l_pac > 50000000? 2 : 3; // set the algorithm for generating BWT
{
strcpy(str, prefix); strcat(str, ".pac");
strcpy(str2, prefix); strcat(str2, ".rpac");
......
......@@ -309,10 +309,6 @@ typedef struct {
bsw2seq1_t *seq;
} bsw2seq_t;
#ifdef HAVE_PTHREAD
static pthread_mutex_t g_dbwtsw_lock = PTHREAD_MUTEX_INITIALIZER;
#endif
static int fix_cigar(const char *qname, const bntseq_t *bns, bsw2hit_t *p, int n_cigar, uint32_t *cigar)
{
// FIXME: this routine does not work if the query bridge three reference sequences
......@@ -469,15 +465,7 @@ static void bsw2_aln_core(int tid, bsw2seq_t *_seq, const bsw2opt_t *_opt, const
l = p->l;
#ifdef HAVE_PTHREAD
if (_opt->n_threads > 1) {
pthread_mutex_lock(&g_dbwtsw_lock);
if (p->tid < 0) p->tid = tid;
else if (p->tid != tid) {
pthread_mutex_unlock(&g_dbwtsw_lock);
continue;
} // in pinciple else should not happen
pthread_mutex_unlock(&g_dbwtsw_lock);
}
if (x % _opt->n_threads != tid) continue;
#endif
// set opt->t
......
#include <stdio.h>
#include <string.h>
#include "main.h"
#include "utils.h"
#ifndef PACKAGE_VERSION
#define PACKAGE_VERSION "0.5.9-r16"
#define PACKAGE_VERSION "0.5.9-r26-dev"
#endif
static int usage()
......@@ -59,5 +60,7 @@ int main(int argc, char *argv[])
fprintf(stderr, "[main] unrecognized command '%s'\n", argv[1]);
return 1;
}
err_fflush(stdout);
err_fclose(stdout);
return 0;
}
......@@ -30,6 +30,7 @@
#include <stdlib.h>
#include <string.h>
#include <zlib.h>
#include <errno.h>
#include "utils.h"
FILE *err_xopen_core(const char *func, const char *fn, const char *mode)
......@@ -80,3 +81,68 @@ void err_fatal_simple_core(const char *func, const char *msg)
fprintf(stderr, "[%s] %s Abort!\n", func, msg);
abort();
}
size_t err_fwrite(const void *ptr, size_t size, size_t nmemb, FILE *stream)
{
size_t ret = fwrite(ptr, size, nmemb, stream);
if (ret != nmemb)
{
err_fatal_simple_core("fwrite", strerror(errno));
}
return ret;
}
int err_printf(const char *format, ...)
{
va_list arg;
int done;
va_start(arg, format);
done = vfprintf(stdout, format, arg);
int saveErrno = errno;
va_end(arg);
if (done < 0)
{
err_fatal_simple_core("vfprintf(stdout)", strerror(saveErrno));
}
return done;
}
int err_fprintf(FILE *stream, const char *format, ...)
{
va_list arg;
int done;
va_start(arg, format);
done = vfprintf(stream, format, arg);
int saveErrno = errno;
va_end(arg);
if (done < 0)
{
err_fatal_simple_core("vfprintf", strerror(saveErrno));
}
return done;
}
int err_fflush(FILE *stream)
{
int ret = fflush(stream);
if (ret != 0)
{
err_fatal_simple_core("fflush", strerror(errno));
}
return ret;
}
int err_fclose(FILE *stream)
{
int ret = fclose(stream);
if (ret != 0)
{
err_fatal_simple_core("fclose", strerror(errno));
}
return ret;
}
......@@ -31,6 +31,15 @@
#include <stdio.h>
#include <zlib.h>
#ifdef __GNUC__
// Tell GCC to validate printf format string and args
#define ATTRIBUTE(list) __attribute__ (list)
#else
#define ATTRIBUTE(list)
#endif
#define err_fatal_simple(msg) err_fatal_simple_core(__func__, msg)
#define xopen(fn, mode) err_xopen_core(__func__, fn, mode)
#define xreopen(fn, mode, fp) err_xreopen_core(__func__, fn, mode, fp)
......@@ -46,6 +55,13 @@ extern "C" {
FILE *err_xopen_core(const char *func, const char *fn, const char *mode);
FILE *err_xreopen_core(const char *func, const char *fn, const char *mode, FILE *fp);
gzFile err_xzopen_core(const char *func, const char *fn, const char *mode);
size_t err_fwrite(const void *ptr, size_t size, size_t nmemb, FILE *stream);
int err_fprintf(FILE *stream, const char *format, ...)
ATTRIBUTE((format(printf, 2, 3)));
int err_printf(const char *format, ...)
ATTRIBUTE((format(printf, 1, 2)));
int err_fflush(FILE *stream);
int err_fclose(FILE *stream);
#ifdef __cplusplus
}
......
#!/usr/bin/perl -w
use strict;
use warnings;
while (<>) {
if (/\tXA:Z:(\S+)/) {
my $l = $1;
print;
my @t = split("\t");
while ($l =~ /([^,;]+),([-+]\d+),([^,]+),(\d+);/g) {
my $mchr = ($t[6] eq $1)? '=' : $t[6]; # FIXME: TLEN/ISIZE is not calculated!
my $seq = $t[9];
my $phred = $t[10];
# if alternative alignment has other orientation than primary,
# then print the reverse (complement) of sequence and phred string
if ((($t[1]&0x10)>0) xor ($2<0)) {
$seq = reverse $seq;
$seq =~ tr/ACGTacgt/TGCAtgca/;
$phred = reverse $phred;
}
print(join("\t", $t[0], 0x100|($t[1]&0x6e9)|($2<0?0x10:0), $1, abs($2), 0, $3, @t[6..7], 0, $seq, $phred, "NM:i:$4"), "\n");
}
} else { print; }
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment