両方のファイルのデータをマージする際に問題があります。これは染色体、位置、参照、および代替対立遺伝子を含む遺伝データであり、どちらも参照および代替対立遺伝子と一緒に4つの列をすべて一致させてファイルをマージする必要があります。したがって、ファイルの$ 1、$ 2、$ 4、および$ 5または$ 1、$ 2、$ 5、および$ 4の列とデータファイルの$ 1、$ 5、$ 6、および$ 7の列が正確に一致するものを見つける必要があります。データファイル内の正確な順序を維持することは非常に重要です。並べ替えることはできません(残念ながら結合を使用しません。これは、このタイプの質問の別のインスタンスで見つかった提案の答えです)。
私はawkを使って数千行のサンプルファイルで作業できるコードを得ましたが、私の大規模なデータセットに合わせて拡張されません(照会ファイルには3億行以上、データファイルには3000万行あります) - おそらくコードこれは、照会のために2つの巨大な配列のメモリを予約する必要があるためです。拡張可能なコード(?Perl?)の提案を送ってくれてありがとう!
検索ファイルの形式は次のとおりです。
1 10150 rs371194064 C T
1 10165 rs796884232 A AC
1 10177 rs367896724 A AC
1 10177 rs201752861 A C
1 10180 rs201694901 T C
1 10199 rs905327004 A T
1 10231 rs200279319 C A
1 10234 rs145599635 C T
1 10235 rs540431307 T TA
1 10235 rs1035249121 T A
1 10235 rs1035249121 T C
1 10241 rs960927773 T C
1 10247 rs796996180 T C
1 10248 rs148908337 A T
1 10249 rs774211241 AAC A
私のデータファイルの形式は次のとおりです。
1 chr1 chr1:10177 1:10177_A_AC 10177 A AC
1 chr1 chr1:10235 1:10235_T_TA 10235 T TA
1 chr1 chr1:10352 1:10352_T_TA 10352 T TA
1 chr1 chr1:10505 1:10505_A_T 10505 A T
1 chr1 chr1:10506 1:10506_C_G 10506 C G
1 chr1 chr1:10511 1:10511_G_A 10511 G A
1 chr1 chr1:10539 1:10539_C_A 10539 C A
1 chr1 chr1:10542 1:10542_C_T 10542 C T
1 chr1 chr1:10579 1:10579_C_A 10579 C A
出力は次のようになります。
1 rs367896724 1:10177_A_AC 10177 A AC A AC
1 rs540431307 1:10235_T_TA 10235 T TA T TA
1 chr1:10352 1:10352_T_TA 10352 T TA T TA
1 chr1:10505 1:10505_A_T 10505 A T A T
1 chr1:10506 1:10506_C_G 10506 C G C G
1 chr1:10511 1:10511_G_A 10511 G A G A
1 chr1:10539 1:10539_C_A 10539 C A C A
1 chr1:10542 1:10542_C_T 10542 C T C T
1 chr1:10579 1:10579_C_A 10579 C A C A
サンプルファイルの操作に成功したawkコードは次のとおりです。
awk 'BEGIN {OFS = "\t"} NR==FNR { #lookup file (323 million rows) key = $1 "," $2 "," $4 "," $5 present[key] = 1 ID[key] = $3 Ref[key] = $4 Alt[key] = $5 key1 = $1 "," $2 "," $4 "," $5 present1[key1] = 1 ID1[key1] = $3 Ref1[key1] = $4 Alt1[key1] = $5 next } { # my data file (3 million rows) key = $1 "," $5 "," $6 "," $7 key1 = $1 "," $5 "," $7 "," $6 if (present[key]) print $1, ID[key], $4, $5, $6, $7, Ref[key], Alt[key]; else if (present1[key1]) print $1, ID1[key1], $4, $5, $6, $7, Ref1[key1], Alt1[key1]; else print $1, $3, $4, $5, $6, $7, $6, $7 }' $lookupfile $mydatafile > $outputfile
答え1
問題は、ファイルをメモリに保存するのではなく、データファイルの各行の検索テーブルを取得することです。あなたのコードには表示されませんが、後ろに3'000'000 x 323'000'000/2 =ほぼ1000兆の文字列比較を実行して、メモリバスから数千テラバイトを移動します。 200GBit / sの高速メモリを使用しても数時間かかります。
したがって、問題の鍵は、ルックアップテーブルを保存する方法です。実行時間を指数関数的に短縮するためにバイナリツリーを使用することをお勧めします。他の言語perl
を使用してC
これを行うことができますが、この時点ではトピックから逸脱しました。
unixコマンドツールセットはこの問題を解決するのに役立ちません。
答え2
私の仮説が正しい場合、両方のファイルは染色体、塩基対位置、rs番号(照会テーブルのみ)、最後に対立遺伝子でソートされます。少なくとも表示されている部分はこのパターンに従います。この場合、照会テーブル全体をメモリに保持する必要はありません。代わりに、無視できるメモリ要件で各ファイルを一度だけ検索するだけです。
データファイルの各タグを順番に繰り返して、一致が見つかるか、候補位置を超えて一致が決定しなくなるまで検索ファイルから検索します。一致するものが見つかると、対応するrs番号がルックアップテーブルから抽出され、それ以外の場合はデータテーブルの現在のchr:bpの組み合わせのみが使用されます。
以下のスクリプトを使用して目的の結果を得ました。スクリプトを保存し、次のように使用します。
gawk -f scriptname datafile lookuptable outputfile
いくつかの小さな追加:処理されたデータ量に関する最小限のフィードバックを得るには、「#」と「。」を使用します。これは、それぞれデータテーブルとルックアップテーブルの10,000行ごとに出力されます。
#!/usr/bin/gawk -f
BEGIN {
OFS = "\t"
step = 10000
while (1==1) {
if ((getline indata < ARGV[1]) < 1)
break
if (!(na++ % step))
printf "\n#"
split(indata,a)
allequal = 0
while (1==1) {
if (!overrun) {
if ((getline inlookup < ARGV[2]) < 1)
break
if (!(nb++ % step))
printf "."
} else {
overrun=0
}
split(inlookup,b)
if (b[1]>a[1] || b[2]>a[5]) {
overrun=1
break
}
if (a[1]==b[1] && a[5]==b[2] && ((a[6]==b[4] && a[7]==b[5]) || (a[7]==b[4] && a[6]==b[5]))) {
allequal=1
break
}
}
if (allequal) {
print a[1],b[3],a[4],a[5],a[6],a[7],b[4],b[5] > ARGV[3]
} else {
print a[1],a[3],a[4],a[5],a[6],a[7],a[6],a[7] > ARGV[3]
}
}
}