2

現在、ファイルに遺伝子のリストがあります。各行には、その情報を含む染色体があります。このようなエントリは次のように表示されます。

NM_198212 chr7 + 115926679 115935830 115927071 11593344 2 115926679、'115933260'、 115927221、'115935830'、

染色体の配列は、塩基 115926679 から始まり、塩基115935830まで (ただし含まない)続きます。

スプライシングされた配列が必要な場合は、エクソンを使用します。最初は 115926679 から 155927221に拡張され、2 番目は '115933260' から '115935830' に拡張されます。

ただし、次のような補完的なシーケンスで問題が発生しました。

NM_001005286 chr1 - 245941755 245942680 245941755 245942680 1 245941755、「245942680」

列 3 は「-」であるため、これらの座標はアンチセンス鎖 (鎖の補体) を参照しています。最初の塩基 (太字) は、センス鎖の最後の塩基 (イタリック体) と一致します。ファイルにはセンス スタンドしかないため、アンチセンス鎖の座標をセンス鎖に変換し、正しい配列を選択してから逆補完する必要があります。

とは言っても、私はプログラミングを始めてまだ半年ほどで、どうやって始めればいいのかわかりません。

私は正規表現を書きました：

'(NM_\d+)\s+(chr\d+)([(\+)|(-)])\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+),(\d+),s+(\d+),(\d+),'

しかし、この機能を開始する方法がわかりません...誰かがこれを始めるのを手伝ってくれて、おそらくこれを行う方法を教えてくれたら、とても感謝しています。

OK: これが 25 番染色体だとします:

AAAAAAAAAACCCCCCCCCCTTTTTTTTTTGGGGGGGGGG

（各キャラクターは10人です）。

今:スプライスされていない遺伝子を探している場合: chr25 + 10 20

次に、遺伝子は 10 番目の位置 (0 から開始) から始まり、20 番目の位置まで上昇しますが、20 番目の位置は含まれません。したがって、次のようになります。

CCCCCCCCCC

かんたんだよ。Python の文字列スライスと非常によく一致します。

私があなたに与えると、より混乱します：

chr25 - 10 20

あなたが持っているのはポジティブストランドです。しかし、この遺伝子はマイナス（相補）鎖にあります。染色体が二本鎖としてどのように見えるかを思い出してください。

AAAAAAAAAACCCCCCCCCCTTTTTTTTTTGGGGGGGGGGGGGGGGGGAAAAAAAAAACCCCCCCCCC
_

一番下の鎖の遺伝子を探しています。つまり、右から 0 から数えます。上のストランドには左から、下のストランドには右から番号を付けます。ここで欲しいのは AAAAAAAAAA です。

キャッチは、私があなたにトップストランドだけを与えているということです. 私はあなたに一番下のストランドを与えていません。(一番上のストランドから自分自身を生成することもできますが、それがどれほど大きいかを考えると、そうしないことをお勧めします。)

したがって、座標を変換する必要があります。下の鎖では、塩基 0 (一番右の C) が上の鎖の塩基 39 と反対になっています。底 1 は底 38 に対して、底 2 はケース 37 に対してです。

つまり、下の鎖の塩基 10-20 を見つけたい場合は、上の鎖の塩基 20-29 を見ることができます (そして、それを逆補数します)。

下部ストランドの座標を下部ストランドの同等の座標に変換する方法を理解する必要があります。はい: とても紛らわしいです

私はウェロニカの元の答えを試しました：

fields = line.split(' \t')
geneID, chr, strand = fields[:2]
start = int(fields[3])
end = int(fields[4])
if strand == '-':
    start,end = -(start + 1), -(end + 1) # this part was changed from original answer.

これは正しい軌道に乗っていますが、十分ではありません。これは 10 と 20 を取り、それを 20 と 10 に変えます。

そして、これを行うことで文字列を逆補完できることを知っています：

r = s[::-1]
bc = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
l = list(r)
o = [bc[base] for base in l]
     return ''.join(o)

編集しました！これは正しいように見えますか?!

fp2 = open('chr22.fa', 'r')
fp = open('chr22.fa', 'r')
for line in fp2:
    newstring = ''
    z = line.strip()
    newstring += z
for line in fp:
    fields = line.split('\t')
    gene_ID, chr, strand = fields[:2]
    start = int(fields[3])
    end = int(fields[4])
    bc = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'a': 't', 't': 'a', 'c':'g', 'g':'c', 'N':'N', 'n':'n'}
    l = list(newstring)        
    if strand == '+':
        geneseq = ''.join([bc[base] for base in l[start:end]]) 
    if strand == '-':
        newstart, newend = -(start + 1), -(end + 1)
        genseq = ''.join([bc[base] for base in l[newstart:newend]])
4

4 に答える 4

0

（特にファイルが大きいので）ファイルバッファから直接読み書きする方がはるかに簡単だと思います。

ヘッダーファイルをすでに解析したとしましょう。解析した行は次のようになります。

line = "NM_001005286 chr1 - 245941755 245942680 245941755 245942680 1"

次に、開始/終了位置が（アンチセンス座標で）何であるかを決定します。

name, chromosome, direction, start, end = line[:5]

次に、次の手順を実行します。

#Open up the file `chr1.txt`.
f = open("chr1.txt", "r")

#Determine the read length.
read_len = end - start

#Seek to the appropriate position (notice the second argument, 2 -- this means
#seek from the end of the file)
f.seek(-end, 2)

#Read the data
sense_seq = f.read(read_len)

その後は、配列をアンチセンスに変換するだけです。

簡単 ：）

于 2012-04-23T01:15:21.523 に答える
0

質問に対する私のコメントで述べたように、これは非常に奇妙なファイル形式のように思われるため、最初は混乱しました。

注: これが標準的な生物学のファイル形式の 1 つである場合は、 Biopythonなどを使用して解析した方がよいでしょう。

独自の解析を行いたい場合、正規表現は依然として間違った方法のように思えます。単純なスペース/タブ区切りのファイルでは読みにくく、不要です。染色体 fasta ファイルを解析し、すべての染色体のシーケンスをchrom_sequencesname:seq 辞書として持っていること、およびreverse_complement関数があることを前提としています (これらはどちらも手作業で簡単に実装できますが、おそらく biopython を使用したほうがよいでしょう)。 .

fields = line.split(' ')  # or '\t' instead of ' ' if the file is tab-separated
gene_ID, chr, strand = fields[:2]
start, end = [int(x) for x in fields[3:5])
this_chromosome_seq = chrom_sequences[chr]
# if strand is +, just take the sequence based on the start-end position
if strand == '+':
  # be careful about off-by-one errors here - some formats give you a 1-based position, 
  #  other formats make it 0-based, and they can also be end-inclusive or end-exclusive
  gene_sequence = this_chromosome_seq[start:end]
# if your coordinates really are given as antisense strand coordinates when strand is -,
#  you just need to subtract them from the chromosome length to get sense-strand coordinates,
#  (switching start and end so they're still in smaller-to-larger order),
#  and then reverse-complement the resulting sequence. 
if strand == '-':
  chrom_length = len(this_chromosome_seq)
  # again, be careful about off-by-one errors here!
  new_start,new_end = chrom_length-end, chrom_length-start
  gene_sequence = reverse_complement(this_chromosome_seq[new_start:new_end])

元の答え、実際には尋ねられたことをしませんでした：

開始位置と終了位置を取得するだけの場合は、次のようにします。

fields = line.split(' ')  # or '\t' instead of ' ' if the file is tab-separated
gene_ID, chr, strand = fields[:2]
start = int(fields[3])
end = int(fields[4])
if strand == '-':
  start,end = end,start

次に、 fasta ファイルを解析して、これらの始点と終点の座標のシーケンスを実際に取得し、それを逆に補完する必要がありstrand=='-'ます. 繰り返しますが、Biopython はそのほとんどを実行できると思います。

別のメモ - Biostarは、バイオインフォマティクスに特化した StackExchange タイプの優れたサイトです。より良い回答が得られる可能性があります。

于 2012-04-21T22:42:04.927 に答える
0

負の数で文字列をスライスすると、Python は末尾から逆方向にカウントします。

complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
s = "AAAAAAAAAACCCCCCCCCCTTTTTTTTTTGGGGGGGGGG"
"".join(complement[c] for c in s[-20:-10])

編集：

はい、あなたの編集は私にはほぼ正しいように見えます。私はフェンスポストのエラーをチェックするのがとても苦手ですが、とにかく私よりもあなたの方がそれらを見るのに適しています!

あなたのコードをより Pythonic になるように書き直しましたが、実質的な変更はありません。

bc = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'N':'N'}

f = open('chr22.fa')
l = ''.join(line.strip() for line in f)
f.seek(0)

for line in f:
    fields = line.split('\t')
    gene_ID, chr, strand = fields[:2]

    start = int(fields[3])
    end = int(fields[4])

    if strand == '-':
        start, end = -(start + 1), -(end + 1)

    geneseq = ''.join(bc[base.upper()] for base in l[start:end])
于 2012-04-22T16:16:20.543 に答える
0

ドメインの問題はわかりませんでしたが、1 つの正規表現に詰め込みすぎているようです。次のように (疑似コードで)、より単純なサブ問題に分割してみてください。

if third column is a '+'
    parseRegularSequence()
else
    parseComplementarySequence()
于 2012-04-21T21:51:22.620 に答える