Mercurial > repos > ryanmorin > nextgen_variant_identification
comparison SNV/SNVMix2_source/SNVMix2-v0.12.1-rc1/samtools-0.1.6/bam_mate.c @ 0:74f5ea818cea
Uploaded
| author | ryanmorin |
|---|---|
| date | Wed, 12 Oct 2011 19:50:38 -0400 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:74f5ea818cea |
|---|---|
| 1 #include <stdlib.h> | |
| 2 #include <string.h> | |
| 3 #include "bam.h" | |
| 4 | |
| 5 // currently, this function ONLY works if each read has one hit | |
| 6 void bam_mating_core(bamFile in, bamFile out) | |
| 7 { | |
| 8 bam_header_t *header; | |
| 9 bam1_t *b[2]; | |
| 10 int curr, has_prev; | |
| 11 | |
| 12 header = bam_header_read(in); | |
| 13 bam_header_write(out, header); | |
| 14 | |
| 15 b[0] = bam_init1(); | |
| 16 b[1] = bam_init1(); | |
| 17 curr = 0; has_prev = 0; | |
| 18 while (bam_read1(in, b[curr]) >= 0) { | |
| 19 bam1_t *cur = b[curr], *pre = b[1-curr]; | |
| 20 if (has_prev) { | |
| 21 if (strcmp(bam1_qname(cur), bam1_qname(pre)) == 0) { // identical pair name | |
| 22 cur->core.mtid = pre->core.tid; cur->core.mpos = pre->core.pos; | |
| 23 pre->core.mtid = cur->core.tid; pre->core.mpos = cur->core.pos; | |
| 24 if (pre->core.tid == cur->core.tid && !(cur->core.flag&(BAM_FUNMAP|BAM_FMUNMAP)) | |
| 25 && !(pre->core.flag&(BAM_FUNMAP|BAM_FMUNMAP))) | |
| 26 { | |
| 27 uint32_t cur5, pre5; | |
| 28 cur5 = (cur->core.flag&BAM_FREVERSE)? bam_calend(&cur->core, bam1_cigar(cur)) : cur->core.pos; | |
| 29 pre5 = (pre->core.flag&BAM_FREVERSE)? bam_calend(&pre->core, bam1_cigar(pre)) : pre->core.pos; | |
| 30 cur->core.isize = pre5 - cur5; pre->core.isize = cur5 - pre5; | |
| 31 } else cur->core.isize = pre->core.isize = 0; | |
| 32 if (pre->core.flag&BAM_FREVERSE) cur->core.flag |= BAM_FMREVERSE; | |
| 33 else cur->core.flag &= ~BAM_FMREVERSE; | |
| 34 if (cur->core.flag&BAM_FREVERSE) pre->core.flag |= BAM_FMREVERSE; | |
| 35 else pre->core.flag &= ~BAM_FMREVERSE; | |
| 36 if (cur->core.flag & BAM_FUNMAP) { pre->core.flag |= BAM_FMUNMAP; pre->core.flag &= ~BAM_FPROPER_PAIR; } | |
| 37 if (pre->core.flag & BAM_FUNMAP) { cur->core.flag |= BAM_FMUNMAP; cur->core.flag &= ~BAM_FPROPER_PAIR; } | |
| 38 bam_write1(out, pre); | |
| 39 bam_write1(out, cur); | |
| 40 has_prev = 0; | |
| 41 } else { // unpaired or singleton | |
| 42 pre->core.mtid = -1; pre->core.mpos = -1; pre->core.isize = 0; | |
| 43 if (pre->core.flag & BAM_FPAIRED) { | |
| 44 pre->core.flag |= BAM_FMUNMAP; | |
| 45 pre->core.flag &= ~BAM_FMREVERSE & ~BAM_FPROPER_PAIR; | |
| 46 } | |
| 47 bam_write1(out, pre); | |
| 48 } | |
| 49 } else has_prev = 1; | |
| 50 curr = 1 - curr; | |
| 51 } | |
| 52 if (has_prev) bam_write1(out, b[1-curr]); | |
| 53 bam_header_destroy(header); | |
| 54 bam_destroy1(b[0]); | |
| 55 bam_destroy1(b[1]); | |
| 56 } | |
| 57 | |
| 58 int bam_mating(int argc, char *argv[]) | |
| 59 { | |
| 60 bamFile in, out; | |
| 61 if (argc < 3) { | |
| 62 fprintf(stderr, "samtools fixmate <in.nameSrt.bam> <out.nameSrt.bam>\n"); | |
| 63 return 1; | |
| 64 } | |
| 65 in = (strcmp(argv[1], "-") == 0)? bam_dopen(fileno(stdin), "r") : bam_open(argv[1], "r"); | |
| 66 out = (strcmp(argv[2], "-") == 0)? bam_dopen(fileno(stdout), "w") : bam_open(argv[2], "w"); | |
| 67 bam_mating_core(in, out); | |
| 68 bam_close(in); bam_close(out); | |
| 69 return 0; | |
| 70 } |