問題タブ [sequence-alignment]

2 つのシーケンス間のローカル アラインメントを見つけるコードを書いています。これは、私が取り組んできた最小限の実用的な例です。

次のコードは、出力で実行されます

しかし、次のように、5 つのアラインメントの横にある「CC」を見つけて、「6」のスコアが可能であるはずです。

何が起こっているかについてのアイデアはありますか？

概要：

整列アルゴリズムの詳細から始めて、最後に質問をします。整列アルゴリズムについて知っている場合は、最初に渡します。

次のような 2 つの文字列があるとします。

Smith-WatermanまたはNeedleman-Wunschのようないくつかのアルゴリズムがあり、この 2 つのシーケンスを整列させて行列を作成します。次のセクションで結果を見てください。

質問：

私の質問は単純ですが、答えは見た目ほど簡単ではないかもしれません。のように文字のグループを単一のものとして使用したい: [A0][C0][A1][B1]. しかし、これらのアルゴリズムでは、個々の文字を使用する必要があります。どうすればそれを達成できますか?

PS 次のシーケンスがあるとします#read #write #add #write。次に、これを次のようなものに変換します: #read to A .... #write to B .... #add to C. 次に、シーケンスは to: になりますABCB。でも、で始まる言葉がたくさんあります#。また、それらすべてを変換するには ASCII テーブルだけでは不十分です。それから、もっとキャラが必要です。唯一の方法は、[A0] ... [Z9]単語ごとに次のようなものを使用することです。または数字を使用します。

PS: Smith-Waterman のサンプル コードがこのリンクにあります。

PS：そのようなものが欲しい別の投稿がありますが、私が欲しいのは異なります。この質問では、a で始まり で[終わる文字のグループがあります]。eeまた、 is equal toのようなセマンティックを使用する必要はありませんi。

解決するアルゴリズム、または少なくとも次の問題の適切な名前を探しています。

ビット文字列のセットBがあります。アルゴリズムは、次のような最小の (「ビット セットが最も少ない」と定義される) ビット文字列Sを見つける必要があります。

Bのすべてのb に対して、となるシフトN ( ℤ</a> 内) が存在します。(S << N) & b == b

これは、ある種の複数の配列アラインメントの問題のように見えます。

B ( i = 1 .. | B |)の各b _iに対して各N _iを見つけることができれば、 Sはすべてのビットごとの論理和 ( b _i >> N _i ) のように見えます。

私の直感では、最初のステップは、Bに別のビット文字列cが存在し、そのようなシフトMが存在するBからすべてのbを削除することです。次は何ですか？b & (c << M) == b

12 種の多数のアラインメント ブロックを含む 600 Mb の FASTA ファイルがあり、それらを対応するアラインメントを持つ 1 つのブロックを含む小さな FASTA ファイルに分割したい

次のようなsedスクリプトがあります。

これは小規模では機能しますが、600Mb のような大きなファイルの場合、約 2 日かかりすぎます。これは私が実行しているコンピューターの問題ではないと思います。

これをより速くする方法を知っている人はいますか？

入力 Fasta ファイルは次のようになります。

dm3.chr3R(-):17092630-17092781|sequence_index=0|block_index=4|species=dm3|dm3_4_0 GGCGGAGATCAAGAATCGCGTCGGGCCGCCGTCCAGCGCCACTGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAACACCAAATCCGGTTGCTGCTGAATAA-CTCTGATTGTGAATCattattttattatacaatta droGri2.scaffold_15074(-):2610183-2610334|sequence_index=0|block_index=4|species=droGri2| droGri2_4_0 GGCGGAGATCAAGAATCGTGTTGGGCCGCCGTCGAGCGCCACCGATAACGCTAGCAAAGTGAAAATCGATCAGGGACGCCCAGTGGAAAACAATAGATCTGGTTGCTGCTAAATAA-CTCTGATTGTGAATCATTATTTTATTATACAATTa droMoj3.scaffold_6540(+):33866311-33866462|sequence_index=0|block_index=4|species=droMoj3|droMoj3_4_0 TGCCGAGATTAAGAATCGTGTCGGTCCGCCGTCCAGCGCAACCGACAATGCAAGCAAAGTGAAAATCGATCAGGGACGTCCAGTGGAGAACACCAGATCTGGTTGCTGCTGAATAA-CTCTGATTGTGAATCATTATTTTATTatacaatta droVir3.scaffold_12822(+):1248119-1248270|sequence_index=0|block_index=4|species=droVir3|droVir3_4_0 GGCCGAGATTAAGAATCGCGTCGGGCCGCCGTCCAGCGCCACCGATAATGCTAGCAAAGTGAAAATCGATCAGGGTCGTCCAGTGGAGAACACCAAATCTGGTTGCTGCTGAATAA-CTCTGATTGTGAATCattattttattatacaatta droWil1.scaffold_181130(-):16071336-16071488|sequence_index=0|block_index=4|species=droWil1|droWil1_4_0 GGCCGAGATTAAGAATCGTGTTGGGCCGCCGTCCAGCGCCACTGATAATGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAATACCAAATCCGGTTGCTGCTGAATAAACTCTGATTGTGAATCATTATTTTATTATACAATTA droPer1.super_19(- ):1310088-1310239|sequence_index=0|block_index=4|species=droPer1|droPer1_4_0 GGCTGAGATCAAGAATCGCGTCGGACCGCCGTCCAGCGCCACCGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAAACCCAATTCTGGTTGCTGCTGAATAA-CTCTGATTGTGAATCattattttatatacaatta dp4(-2):chr4.chr.5593491-5593642|sequence_index=0|block_index=4|species=dp4|dp4_4_0 GGCTGAGATCAAGAATCGCGTCGGACCGCCGTCCAGCGCCACCGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAAGCCCAATTCTGGTTGCTGCTGAATAA-CTCTGATTGTGAATCattattttattatacaatta droAna3.scaffold_13340(-):3754154-3754305|sequence_index=0|block_index=4|species=droAna3|droAna3_4_0 GGCCGAGATCAAGAATCGCGTCGGGCCACCGTCCAGCGCCACCGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAACACCAGATCCGGTTGCTGCTGAATAA-CTCTGATTGTGAATCattattttattataaaatta droEre2.scaffold_4770 (+):4567591-4567742|sequence_index=0|block_index=4|species=droEre2|droEre2_4_05883047-5883198|sequence_index=0|block_index=4|species=droYak2|droYak2_4_0 GGCCGAGATCAAGAATCGCGTCGGGCCGCCATCCAGCGCCACCGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAACACCAAATCCGGTTGCTGCTGAATAA-CTCTGATTGTGAATCattattttattatacaatta droSec1.super_38(+):36432-36583|sequence_index=0|block_index=4|species=droSec1|droSec1_4_0 GGCGGAGATCAAGAATCGCGTCGGTCCGCCGTCCAGCGCCACTGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAACACCAAATCCGGTTGCTGCTGAATAA-CTCTGATTGTGAATCattattttattatacaatta droSim1.chr3R (+):4366350-4366501|sequence_index=0|block_index=4|species=droSim1|droSim1_4_04366350-4366501|sequence_index=0|block_index=4|種=droSim1|droSim1_4_04366350-4366501|sequence_index=0|block_index=4|種=droSim1|droSim1_4_0super_38(+):36432-36583|sequence_index=0|block_index=4|species=droSec1|droSec1_4_0 GGCGGAGATCAAGAATCGCGTCGGTCCGCCGTCCAGCGCCACTGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAACACCAAATCCGGTTGCTGCTGAATAA-CTCTGATTGTGAATCattattttattatacaatta droSim1.chr3R(+):4366350-4366501|sequence_index=0|block_index=4|species=droSim1|droSim1_4_0 GGCGGAGATCAAGAATCGCGTCGGGCCGCCGTCCAGCGCCACTGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAACACCAAATCCGGTTGCTGCTGAATAA -CTCTGATTGTGAATCattattttattatacaattasuper_38(+):36432-36583|sequence_index=0|block_index=4|species=droSec1|droSec1_4_0 GGCGGAGATCAAGAATCGCGTCGGTCCGCCGTCCAGCGCCACTGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAACACCAAATCCGGTTGCTGCTGAATAA-CTCTGATTGTGAATCattattttattatacaatta droSim1.chr3R(+):4366350-4366501|sequence_index=0|block_index=4|species=droSim1|droSim1_4_0 GGCGGAGATCAAGAATCGCGTCGGGCCGCCGTCCAGCGCCACTGACAACGCTAGCAAAGTGAAAATCGATCAAGGACGTCCAGTGGAAAACACCAAATCCGGTTGCTGCTGAATAA -CTCTGATTGTGAATCattattttattatacaatta4366350-4366501|sequence_index=0|block_index=4|種=droSim1|droSim1_4_04366350-4366501|sequence_index=0|block_index=4|種=droSim1|droSim1_4_0

dm3.chr3R(-):17092781-17092867|sequence_index=0|block_index=5|species=dm3|dm3_5_0 GAGTACGCCGCCCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAACGTTGAGCAGGCCTTCATGACGATGGC droSim1.chr3R(+):4366264-4366350|sequence_index=0|block_index=5|species=droSim1|droSim1_5_0 GAGTACGCCGCCCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAACGTTGAGCAGGCCTTTATGACGATGGC droSec1.super_38(+):36346-36432|sequence_index=0|block_index=5|species=droSec1|droSec1_5_0 GAGTACGCCGCCCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAACGTTGAGCAGGCCTTCATGACGATGGC droYak2.chr3R(-):5883198-5883284|sequence_index=0|block_index=5|species=droYak2|droYak2_5_0 GAGTACGCCGCCCAGTTAGGCATTCCATTCCTTGAAACATCGGCCAAGAGCGCCACCAACGTGGAGCAGGCCTTCATGACGATGGC droEre2.scaffold_4770(+):4567505-4567591|sequence_index=0|block_index=5|species=droEre2|droEre2_5_0 GAGTACGCCGCCCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAACGTGGAGCAGGCCTTCATGACGATGGC droAna3.scaffold_13340(+):20375068-20375148|sequence_index=0|block_index=5|species=droAna3|droAna3_5_0 ------GCCGAAAACTTCGACATGCCCTTCTTCGAGGTCTCTTGCAAGTCAAACATCAATATTGAAGATGCGTTTCTTTCCCTGGC dp4 .chr2(-):5593642-5593728|sequence_index=0|block_index=5|species=dp4|dp4_5_0 GAGTATGCAGCTCAGTTAGGCATTCCATTTCTTGAAACTTCGGCCAAGAGCGCCACGAACGTGGAGCAGGCCTTCATGACGATGGC droPer1.super_19(-):1310239-1310325|sequence_index=0|block_index=5|species=droPer1|droPer1_5_0 GAGTATGCAGCTCAGTTAGGCATTCCATTTCTTGAAACTTCGGCCAAGAGCGCCACGAACGTGGAGCAGGCCTTCATGACGATGGC droWil1 .scaffold_181130(-):16071488-16071574|sequence_index=0|block_index=5|species=droWil1|droWil1_5_0 GAATATGCGGCTCAGTTAGGCATTCCATTCCTTGAAACTTCGGCAAAGAGTGCCACCAATGTGGAGCAGGCCTTTATGACGATGGC droVir3.scaffold_12822(+):1248033-1248119|sequence_index=0|block_index=5|species=droVir3|droVir3_5_0 GAGTACGCACATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAACGTGGAGCAGGCATTTATGACGATGGC droMoj3.scaffold_6540(+): 33866225-33866311|sequence_index=0|block_index=5|species=droMoj3|droMoj3_5_0 GAGTATGCACATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAATGTAGAGCAGGCATTCATGACGATGGC droGri2.scaffold_15074(-):2610334-2610420|sequence_index=0|block_index=5|species=droGri2|droGri2_5_0 GAGTACGCAAATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCGAAGAGTGCCACCAATGTGGAACAGGCATTCATGACGATGGC1248033-1248119|sequence_index=0|block_index=5|species=droVir3|droVir3_5_0 GAGTACGCACATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAACGTGGAGCAGGCATTTATGACGATGGC droMoj3.scaffold_6540(+):33866225-33866311|sequence_index=0|block_index=5|species=droMoj3|droMoj3_5_0 GAGTATGCACATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAATGTAGAGCAGGCATTCATGACGATGGC droGri2.scaffold_15074(-): 2610334-2610420|sequence_index=0|block_index=5|species=droGri2|droGri2_5_0 GAGTACGCAAATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCGAAGAGTGCCACCAATGTGGAACAGGCATTCATGACGATGGC1248033-1248119|sequence_index=0|block_index=5|species=droVir3|droVir3_5_0 GAGTACGCACATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAACGTGGAGCAGGCATTTATGACGATGGC droMoj3.scaffold_6540(+):33866225-33866311|sequence_index=0|block_index=5|species=droMoj3|droMoj3_5_0 GAGTATGCACATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAATGTAGAGCAGGCATTCATGACGATGGC droGri2.scaffold_15074(-): 2610334-2610420|sequence_index=0|block_index=5|species=droGri2|droGri2_5_0 GAGTACGCAAATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCGAAGAGTGCCACCAATGTGGAACAGGCATTCATGACGATGGC33866225-33866311|sequence_index=0|block_index=5|species=droMoj3|droMoj3_5_0 GAGTATGCACATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAATGTAGAGCAGGCATTCATGACGATGGC droGri2.scaffold_15074(-):2610334-2610420|sequence_index=0|block_index=5|species=droGri2|droGri2_5_0 GAGTACGCAAATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCGAAGAGTGCCACCAATGTGGAACAGGCATTCATGACGATGGC33866225-33866311|sequence_index=0|block_index=5|species=droMoj3|droMoj3_5_0 GAGTATGCACATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCCAAGAGCGCCACCAATGTAGAGCAGGCATTCATGACGATGGC droGri2.scaffold_15074(-):2610334-2610420|sequence_index=0|block_index=5|species=droGri2|droGri2_5_0 GAGTACGCAAATCAGTTAGGCATTCCATTCCTTGAAACTTCGGCGAAGAGTGCCACCAATGTGGAACAGGCATTCATGACGATGGC

120 塩基の約 13,000 のシーケンスがあり、それらを比較して、保存された領域、それらの間の平均的な相違、または非常に異なる外れ値などを見つけたいと考えています。

問題は、この数のシーケンスでは、私が試したことが実行できないことです。

誰かがこのサイズで似たようなことをしたことがありますか?それを達成する方法のヒントを教えてもらえますか? それとも、私が探すべきいくつかのヒントですか？

私はすでにここや他のフォーラムを見てきましたが、私の質問に対する答えが見つかりませんでした. ターゲット エンリッチメント シーケンシング アプローチ用のベイトを設計し、参照として

A. thaliana を使用して、4 つの異なるゲノムからオーソロガス遺伝子座を MarkerMiner 検索の出力を得たいと考えています。
これらの出力アラインメントは、A. thaliana の注釈付き遺伝子ごとに個別の Fasta ファイルであり、私のデータセットからの配列がアラインされています。
スクリプトを実行して、4 つの入力データセットのうち少なくとも 2 つによってオーソロガスであることがサポートされている遺伝子座を除外しました。

しかし、今、私は困惑しています。
入力データはほとんどが RNAseq であるのに対し、リファレンスにはイントロンも含まれているため、私のアラインメントはギャップがあります。したがって、次のようになります。

AT01G1234567 ATCGATCGATGCGCGCTAGCTGAATCGATCGGATCGCGGTAGCTGGAGCTAGSTCGGATCGC MyData1
CGATGCGCGC-----------CGGATCGCGG---------------CGGATCGC
MyData2
CGCTGCGCGC------------GGATAGCGG--- ------------CGGATCCC

ベイトを効果的に設計するには、整列したすべてのパーツをファイルから抽出する必要があります。これにより、個別のファイルが作成されます。またはファイル内の個別の配置。MyData と Reference シーケンスの間でアラインされている部分で、ギャップのあるすべての部分が除外されています。これらの fasta ファイルは約 1300 個あるため、手動で行うことはできません。Python と Linux コマンド ライン ツールを使用したプログラミングの経験は少しありますが、これを行う方法が完全にわかりません。使用できるツールの種類や、考え出す必要があるアルゴリズムの種類について、ヒントをいただければ幸いです。ありがとうございました。乾杯

Needleman-Wunsch アルゴリズム (通常、ヌクレオチド/タンパク質配列の整列に使用されます) の結果をどのように定量化できるのでしょうか。

一定のスコアリング スキームと、長さが異なる と の 2 つのシーケンスを考えてみましょS1うS2。と のすべての可能なアラインメントを力ずくで計算S1しS2、最高スコアのアラインメントにはスコアがあるとしxます。もちろん、これは Needleman-Wunsch アプローチよりもかなり複雑です。

Needleman-Wunsch アルゴリズムを使用して配列アラインメントを見つける場合、スコアがあると言いますy。

は、2 つのランダム シーケンスおよびrに対して Needleman-Wunsch によって生成されたスコアであると考えてください。R1R2

xと比べてどうyですか？相同性が既知の 2 つの配列yよりも常に大きいですか?r

一般的に、Needleman-Wunsch アルゴリズムを使用して配列アラインメントを大幅に高速化することは理解していますが (力ずくのアプローチと比較して)、それに伴う精度のコスト (もしあれば) については理解していません。元の論文 (Needleman & Wunsch, 1970) を読んでみましたが、まだこの疑問が残っています。