htslibSAM/BAM ファイルを読み取るためにライブラリを使用していますが、完全に機能します。アラインメントを新しい SAM/BAM ファイルに書き戻すこともできます。
たとえば、次のコードはアラインメントの DNA 配列を出力します。
bam1_t *b = ...;
int i;
for (i = 0; i < b->core.l_qseq; ++i) {
printf("%c", seq_nt16_str[bam_seqi(bam_get_seq(b),i)]);
}
質問: クエリ シーケンスを変更するにはどうすればよいですか? たとえば、最初の文字を「T」に変更しますか? bam_get_seq読み取りのシーケンスを返しますが、bam_set_seq関数はありませんか? 理想的には、次のようなものを探しています:
bam_set_seq(b, 'TTTT') # My new DNA sequence
更新の方法がわかれば、その情報を新しい SAM/BAM ファイルに書き込む方法がわかります。