← Back to branch summary

~ubuntu-branches/ubuntu/saucy/bwa/saucy

« back to all changes in this revision

Viewing changes to bwape.c

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2011-02-08 16:11:12 UTC
  • mfrom: (1.1.7 upstream)
  • Revision ID: james.westby@ubuntu.com-20110208161112-1fvt0mtygl8rhajf
Tags: 0.5.9-1
* New usptream release
  - Upstream moved to GitHub (debian/upstream-metadata.yaml).
  - Checked copyright; one file was removed (debian/copyright).
* Migrated the Debian source package to Git (debian/control).
* Use Debhelper 8 (debian/control, debian/compat).

Show diffs side-by-side

added added

removed removed

Lines of Context:
bwape.c
45
45
void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, int mode, int max_top2);
46
46
bntseq_t *bwa_open_nt(const char *prefix);
47
47
void bwa_print_sam_SQ(const bntseq_t *bns);
 
48
void bwa_print_sam_PG();
48
49
 
49
50
pe_opt_t *bwa_init_pe_opt()
50
51
{
643
644
 
644
645
void bwa_sai2sam_pe_core(const char *prefix, char *const fn_sa[2], char *const fn_fa[2], pe_opt_t *popt)
645
646
{
 
647
        extern bwa_seqio_t *bwa_open_reads(int mode, const char *fn_fa);
646
648
        int i, j, n_seqs, tot_seqs = 0;
647
649
        bwa_seq_t *seqs[2];
648
650
        bwa_seqio_t *ks[2];
649
651
        clock_t t;
650
652
        bntseq_t *bns, *ntbns = 0;
651
653
        FILE *fp_sa[2];
652
 
        gap_opt_t opt;
 
654
        gap_opt_t opt, opt0;
653
655
        khint_t iter;
654
656
        isize_info_t last_ii; // this is for the last batch of reads
655
657
        char str[1024];
662
664
        for (i = 1; i != 256; ++i) g_log_n[i] = (int)(4.343 * log(i) + 0.5);
663
665
        bns = bns_restore(prefix);
664
666
        srand48(bns->seed);
665
 
        for (i = 0; i < 2; ++i) {
666
 
                ks[i] = bwa_seq_open(fn_fa[i]);
667
 
                fp_sa[i] = xopen(fn_sa[i], "r");
668
 
        }
 
667
        fp_sa[0] = xopen(fn_sa[0], "r");
 
668
        fp_sa[1] = xopen(fn_sa[1], "r");
669
669
        g_hash = kh_init(64);
670
670
        last_ii.avg = -1.0;
671
671
 
672
672
        fread(&opt, sizeof(gap_opt_t), 1, fp_sa[0]);
673
 
        fread(&opt, sizeof(gap_opt_t), 1, fp_sa[1]);
 
673
        ks[0] = bwa_open_reads(opt.mode, fn_fa[0]);
 
674
        opt0 = opt;
 
675
        fread(&opt, sizeof(gap_opt_t), 1, fp_sa[1]); // overwritten!
 
676
        ks[1] = bwa_open_reads(opt.mode, fn_fa[1]);
674
677
        if (!(opt.mode & BWA_MODE_COMPREAD)) {
675
678
                popt->type = BWA_PET_SOLID;
676
679
                ntbns = bwa_open_nt(prefix);
688
691
 
689
692
        // core loop
690
693
        bwa_print_sam_SQ(bns);
691
 
        while ((seqs[0] = bwa_read_seq(ks[0], 0x40000, &n_seqs, opt.mode & BWA_MODE_COMPREAD, opt.trim_qual)) != 0) {
 
694
        bwa_print_sam_PG();
 
695
        while ((seqs[0] = bwa_read_seq(ks[0], 0x40000, &n_seqs, opt0.mode, opt0.trim_qual)) != 0) {
692
696
                int cnt_chg;
693
697
                isize_info_t ii;
694
698
                ubyte_t *pacseq;
695
699
 
696
 
                seqs[1] = bwa_read_seq(ks[1], 0x40000, &n_seqs, opt.mode & BWA_MODE_COMPREAD, opt.trim_qual);
 
700
                seqs[1] = bwa_read_seq(ks[1], 0x40000, &n_seqs, opt.mode, opt.trim_qual);
697
701
                tot_seqs += n_seqs;
698
702
                t = clock();
699
703
 
714
718
 
715
719
                fprintf(stderr, "[bwa_sai2sam_pe_core] print alignments... ");
716
720
                for (i = 0; i < n_seqs; ++i) {
717
 
                        bwa_print_sam1(bns, seqs[0] + i, seqs[1] + i, opt.mode, opt.max_top2);
718
 
                        bwa_print_sam1(bns, seqs[1] + i, seqs[0] + i, opt.mode, opt.max_top2);
 
721
                        bwa_seq_t *p[2];
 
722
                        p[0] = seqs[0] + i; p[1] = seqs[1] + i;
 
723
                        if (p[0]->bc[0] || p[1]->bc[0]) {
 
724
                                strcat(p[0]->bc, p[1]->bc);
 
725
                                strcpy(p[1]->bc, p[0]->bc);
 
726
                        }
 
727
                        bwa_print_sam1(bns, p[0], p[1], opt.mode, opt.max_top2);
 
728
                        bwa_print_sam1(bns, p[1], p[0], opt.mode, opt.max_top2);
719
729
                }
720
730
                fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC); t = clock();
721
731
 
742
752
 
743
753
int bwa_sai2sam_pe(int argc, char *argv[])
744
754
{
 
755
        extern char *bwa_rg_line, *bwa_rg_id;
 
756
        extern int bwa_set_rg(const char *s);
745
757
        int c;
746
758
        pe_opt_t *popt;
747
759
        popt = bwa_init_pe_opt();
748
 
        while ((c = getopt(argc, argv, "a:o:sPn:N:c:f:A")) >= 0) {
 
760
        while ((c = getopt(argc, argv, "a:o:sPn:N:c:f:Ar:")) >= 0) {
749
761
                switch (c) {
 
762
                case 'r':
 
763
                        if (bwa_set_rg(optarg) < 0) {
 
764
                                fprintf(stderr, "[%s] malformated @RG line\n", __func__);
 
765
                                return 1;
 
766
                        }
 
767
                        break;
750
768
                case 'a': popt->max_isize = atoi(optarg); break;
751
769
                case 'o': popt->max_occ = atoi(optarg); break;
752
770
                case 's': popt->is_sw = 0; break;
754
772
                case 'n': popt->n_multi = atoi(optarg); break;
755
773
                case 'N': popt->N_multi = atoi(optarg); break;
756
774
                case 'c': popt->ap_prior = atof(optarg); break;
757
 
        case 'f': freopen(optarg, "w", stdout); break;
 
775
                case 'f': xreopen(optarg, "w", stdout); break;
758
776
                case 'A': popt->force_isize = 1; break;
759
777
                default: return 1;
760
778
                }
768
786
                fprintf(stderr, "         -n INT   maximum hits to output for paired reads [%d]\n", popt->n_multi);
769
787
                fprintf(stderr, "         -N INT   maximum hits to output for discordant pairs [%d]\n", popt->N_multi);
770
788
                fprintf(stderr, "         -c FLOAT prior of chimeric rate (lower bound) [%.1le]\n", popt->ap_prior);
771
 
        fprintf(stderr, "         -f FILE sam file to output results to [stdout]\n\n");
 
789
        fprintf(stderr, "         -f FILE  sam file to output results to [stdout]\n");
 
790
                fprintf(stderr, "         -r STR   read group header line such as `@RG\\tID:foo\\tSM:bar' [null]\n");
772
791
                fprintf(stderr, "         -P       preload index into memory (for base-space reads only)\n");
773
792
                fprintf(stderr, "         -s       disable Smith-Waterman for the unmapped mate\n");
774
793
                fprintf(stderr, "         -A       disable insert size estimate (force -s)\n\n");
779
798
                return 1;
780
799
        }
781
800
        bwa_sai2sam_pe_core(argv[optind], argv + optind + 1, argv + optind+3, popt);
 
801
        free(bwa_rg_line); free(bwa_rg_id);
782
802
        free(popt);
783
803
        return 0;
784
804
}