ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/bam2fastx.cpp
(Generate patch)
# Line 92 | Line 92
92    bool ismapped=((b->core.flag & BAM_FUNMAP) == 0);
93    if (ismapped && !all_reads) return;
94    if (mapped_only && !ismapped) return;
95 +
96 +  bool isreversed=((b->core.flag & BAM_FREVERSE) != 0);
97 +
98 +  for(i=0;i<(b->core.l_qseq);i++)
99 +    qseq[i] = bam1_seqi(s,i);
100 +  qseq[i] = 0;
101    
102 <  for(i=0;i<(b->core.l_qseq);i++) {
103 <    int8_t v = bam1_seqi(s,i);
104 <    qseq[i] = bam_nt16_rev_table[v];
102 >  // copied from sam_view.c in samtools.
103 >  static const int8_t seq_comp_table[16] = { 0, 8, 4, 12, 2, 10, 9, 14, 1, 6, 5, 13, 3, 11, 7, 15 };
104 >  if (isreversed) {
105 >    int qlen = b->core.l_qseq;
106 >    for(i=0;i<qlen>>1;i++) {
107 >      int8_t t = seq_comp_table[qseq[qlen - 1 - i]];
108 >      qseq[qlen - 1 - i] = seq_comp_table[qseq[i]];
109 >      qseq[i] = t;
110      }
111 <  qseq[i] = 0;
112 <  
111 >    if(qlen&1) qseq[i] = seq_comp_table[qseq[i]];
112 >  }
113 >
114 >  for(i=0;i<(b->core.l_qseq);i++) {
115 >    qseq[i] = bam_nt16_rev_table[qseq[i]];
116 >  }
117 >  
118    fprintf(fout, "@%s\n%s\n",name, qseq);
119    for(i=0;i<(b->core.l_qseq);i++) {
120      qseq[i]=qual[i]+33;
121      }
122    qseq[i]=0;
123 +
124 +  if(isreversed) {
125 +    int qlen = b->core.l_qseq;
126 +    for(i= 0;i<qlen>>1;i++) {
127 +      int8_t t = qseq[qlen - 1 - i];
128 +      qseq[qlen - 1 - i] = qseq[i];
129 +      qseq[i] = t;
130 +    }
131 +  }
132 +  
133    fprintf(fout, "+\n%s\n",qseq);
134   }
135  
# Line 111 | Line 137
137    char *name  = bam1_qname(b);
138    char *s    = (char*)bam1_seq(b);
139    int i;
140 +
141 +  bool ismapped=((b->core.flag & BAM_FUNMAP) == 0);
142 +  if (ismapped && !all_reads) return;
143 +  if (mapped_only && !ismapped) return;
144 +
145 +  bool isreversed=((b->core.flag & BAM_FREVERSE) != 0);
146 +
147 +  for(i=0;i<(b->core.l_qseq);i++)
148 +    qseq[i] = bam1_seqi(s,i);
149 +  qseq[i] = 0;
150 +  
151 +  // copied from sam_view.c in samtools.
152 +  static const int8_t seq_comp_table[16] = { 0, 8, 4, 12, 2, 10, 9, 14, 1, 6, 5, 13, 3, 11, 7, 15 };
153 +  if (isreversed) {
154 +    int qlen = b->core.l_qseq;
155 +    for(i=0;i<qlen>>1;i++) {
156 +      int8_t t = seq_comp_table[qseq[qlen - 1 - i]];
157 +      qseq[qlen - 1 - i] = seq_comp_table[qseq[i]];
158 +      qseq[i] = t;
159 +    }
160 +    if(qlen&1) qseq[i] = seq_comp_table[qseq[i]];
161 +  }
162 +
163    for(i=0;i<(b->core.l_qseq);i++) {
164 <    int8_t v = bam1_seqi(s,i);
116 <    qseq[i] = bam_nt16_rev_table[v];
164 >    qseq[i] = bam_nt16_rev_table[qseq[i]];
165    }
166    qseq[i] = 0;
167    fprintf(fout, ">%s\n%s\n",name, qseq);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines