さまざまなバリアントで 1 つのタンデム リピートを含むプログラム出力があります。モチーフを (文字列で) 検索し、最大「3」の不一致/挿入/削除ですべてのバリアントを見つけるようにプログラムに指示することは可能ですか?
3 に答える
提供された非常に限られた情報を使用して、これを調べます。
まず、短いフレンドリーな社説:
<editorial>
少なくとも、次のことをお願いします。
- リンクや正確な定義を提供せずに、「モチーフ」、「タンデムリピート」、「塩基対」などのドメイン固有の専門用語は控えてください。
- 目標は何か、これまでに何をしたかを述べます。
- 重要: 入力と目的の出力の明確な例を提供してください。
あなたの質問を理解しようとするために、コメントで20の質問をしなければならないSOの潜在的なヘルパーには役に立ちません! 私は、あなたが何を求めているのかを理解することに、答えるよりも多くの時間を費やしました。
</editorial>
次のプログラムは、長さ 1,000 要素の配列内に長さ 5,428 ペアの 2 文字ペアの文字列を生成します。これらをファイルから読み取る可能性が高いことは承知していますが、これは単なる例です。明らかに、ランダムな文字列を任意のソースからの実際のデータに置き換えます。
'AT','CG','TC','CA','TG','GC','GG'
私が使用したのが正当な塩基対の組み合わせであるかどうかはわかりません。(私は生物学を通して眠りました...)テスト用に正当なランダム文字列を生成したい場合は、map
ブロックのペアを正当なペアに編集し、ペアの数に変更してください。7
ポイントの部分文字列offset
が 3 つ以下の違いである場合、配列要素 (スカラー値) は、ハッシュの値部分の無名配列に格納されます。ハッシュのキー部分は、ほぼ一致する部分文字列です。値は、配列要素ではなく、モチーフに関連付けるファイル名、Perl データ参照、またはその他の関連参照にすることができます。
文字列間の文字ごとの違いを見てきましたがforeach my $j (0..$#a1) { $diffs++ unless ($a1[$j] eq $a2[$j]); }
、問題に適した比較ロジックで行を置き換えることにより、確認する必要がある特定のロジックを配置できます。あなたの文字列でがどのように表現されているかわかりませんmismatches/insertions/deletions
ので、それは読者の演習として残します。おそらくAlgorithm::DiffまたはString:: CPAN の Diff ですか?
このプログラムを変更して、固定オフセットで複数の文字列を検索するのではなく、キーボード入力を使用$target
したり、文字列を最初から最後まで検索したりするのは簡単です。$offset
繰り返しますが、あなたの目標が何であるかは本当に明確ではありませんでした...
use strict; use warnings;
my @bps;
push(@bps,join('',map { ('AT','CG','TC','CA','TG','GC','GG')[rand 7] }
0..5428)) for(1..1_000);
my $len=length($bps[0]);
my $s_count= scalar @bps;
print "$s_count random strings generated $len characters long\n" ;
my $target="CGTCGCACAG";
my $offset=832;
my $nlen=length $target;
my %HoA;
my $diffs=0;
my @a2=split(//, $target);
substr($bps[-1], $offset, $nlen)=$target; #guarantee 1 match
substr($bps[-2], $offset, $nlen)="CATGGCACGG"; #anja example
foreach my $i (0..$#bps) {
my $cand=substr($bps[$i], $offset, $nlen);
my @a1=split(//, $cand);
$diffs=0;
foreach my $j (0..$#a1) { $diffs++ unless ($a1[$j] eq $a2[$j]); }
next if $diffs > 3;
push (@{$HoA{$cand}}, $i);
}
foreach my $hit (keys %HoA) {
my @a1=split(//, $hit);
$diffs=0;
my $ds="";
foreach my $j (0..$#a1) {
if($a1[$j] eq $a2[$j]) {
$ds.=" ";
} else {
$diffs++;
$ds.=$a1[$j];
}
}
print "Target: $target\n",
"Candidate: $hit\n",
"Differences: $ds $diffs differences\n",
"Array element: ";
foreach (@{$HoA{$hit}}) {
print "$_ " ;
}
print "\n\n";
}
出力:
1000 random strings generated 10858 characters long
Target: CGTCGCACAG
Candidate: CGTCGCACAG
Differences: 0 differences
Array element: 999
Target: CGTCGCACAG
Candidate: CGTCGCCGCG
Differences: CGC 3 differences
Array element: 696
Target: CGTCGCACAG
Candidate: CGTCGCCGAT
Differences: CG T 3 differences
Array element: 851
Target: CGTCGCACAG
Candidate: CGTCGCATGG
Differences: TG 2 differences
Array element: 986
Target: CGTCGCACAG
Candidate: CATGGCACGG
Differences: A G G 3 differences
Array element: 998
..several cut out..
Target: CGTCGCACAG
Candidate: CGTCGCTCCA
Differences: T CA 3 differences
Array element: 568 926
BioPerl には、この種のルーチンがあると思います。
いずれにせよ、バイオインフォマティクス スタック エクスチェンジである BioStar でこれを尋ねると、より良い回答が得られる可能性があります。
私が perl を学び始めて最初の 2 年間は、tandyman と呼ばれる非常に効率の悪い (しかし機能的な) タンデム リピート ファインダー (私の以前の仕事の会社の Web サイトで利用可能だった) であると今では考えているものを書きました。数年後に、cottonTandy というファジー バージョンを書きました。今日それを書き直すとしたら、グローバル検索にハッシュを使用し (許容されるミスを考慮して)、ローカル検索にパターン マッチングを利用します。
使用方法の例を次に示します。
#!/usr/bin/perl
use Tandyman;
$sequence = "ATGCATCGTAGCGTTCAGTCGGCATCTATCTGACGTACTCTTACTGCATGAGTCTAGCTGTACTACGTACGAGCTGAGCAGCGTACgTG";
my $tandy = Tandyman->new(\$sequence,'n'); #Can't believe I coded it to take a scalar reference! Prob. fresh out of a cpp class when I wrote it.
$tandy->SetParams(4,2,3,3,4);
#The parameters are, in order:
# repeat unit size
# min number of repeat units to require a hit
# allowed mistakes per unit (an upper bound for "mistake concentration")
# allowed mistakes per window (a lower bound for "mistake concentration")
# number of units in a "window"
while(@repeat_info = $tandy->FindRepeat())
{print(join("\t",@repeat_info),"\n")}
このテストの出力は次のようになります (実行には恐ろしい 11 秒かかります)。
25 32 TCTA 2 0.87 TCTA TCTG
58 72 CGTA 4 0.81 CTGTA CTA CGTA CGA
82 89 CGTA 2 0.87 CGTA CGTG
45 51 TGCA 2 0.87 TGCA TGA
65 72 ACGA 2 0.87 ACGT ACGA
23 29 CTAT 2 0.87 CAT CTAT
36 45 TACT 3 0.83 TACT CT TACT
24 31 ATCT 2 1 ATCT ATCT
51 59 AGCT 2 0.87 AGTCT AGCT
33 39 ACGT 2 0.87 ACGT ACT
62 72 ACGT 3 0.83 ACT ACGT ACGA
80 88 ACGT 2 0.87 AGCGT ACGT
81 88 GCGT 2 0.87 GCGT ACGT
63 70 CTAC 2 0.87 CTAC GTAC
32 38 GTAC 2 0.87 GAC GTAC
60 74 GTAC 4 0.81 GTAC TAC GTAC GAGC
23 30 CATC 2 0.87 CATC TATC
71 82 GAGC 3 0.83 GAGC TGAGC AGC
1 7 ATGC 2 0.87 ATGC ATC
54 60 CTAG 2 0.87 CTAG CTG
15 22 TCAG 2 0.87 TCAG TCGG
70 81 CGAG 3 0.83 CGAG CTGAG CAG
44 50 CATG 2 0.87 CTG CATG
25 32 TCTG 2 0.87 TCTA TCTG
82 89 CGTG 2 0.87 CGTA CGTG
55 73 TACG 5 0.75 TAGCTG TAC TACG TACG AG
69 83 AGCG 4 0.81 ACG AGCTG AGC AGCG
15 22 TCGG 2 0.87 TCAG TCGG
ご覧のとおり、インデルと SNP が許可されます。列は、順番に次のとおりです。
- 開始位置
- 停止位置
- コンセンサス配列
- 見つかったユニットの数
- 1 段階評価の品質指標
- スペースで区切られた繰り返し単位
ジャンク/取るに足らない「繰り返し」を出力するパラメーターを (上記の出力からわかるように) 提供するのは簡単ですが、適切なパラメーターを提供する方法を知っていれば、検索時に設定したものを見つけることができます。
残念ながら、パッケージは公開されていません。それは非常に遅く、原核生物サイズのゲノム検索にも適していないため、わざわざ利用できるようにすることはありませんでした (ただし、個々の遺伝子に対しては実行可能ですが)。初心者のコーディング時代に、「状態」を入力として受け取る機能を追加して、シーケンスのセクションで並列に実行できるようにし始めました。その時点で、私は他のプロジェクトに移っていました。しかし、それがあなたのニーズに合っている場合は、私にメッセージを送ってください.コピーをメールで送信できます.
わずか 1000 行のコードですが、IUPAC 曖昧コード (BDHVRYKMSWN) の許容など、多くの付加機能があります。アミノ酸にも核酸にも効きます。内部反復を除外します (たとえば、TTTT または ATAT を 4nt コンセンサスとして報告しません)。