여러 열과 일치하는 2개의 큰 파일을 병합하고 순서를 유지합니다(일치하는 값과 일치하지 않는 값 인쇄) - awk에서 확장됨

여러 열과 일치하는 2개의 큰 파일을 병합하고 순서를 유지합니다(일치하는 값과 일치하지 않는 값 인쇄) - awk에서 확장됨

두 파일의 데이터를 병합하는 데 문제가 있습니다. 이는 염색체, 위치, 참조 및 대체 대립 유전자가 포함된 유전 데이터이며 어느 쪽이든 참조 및 대체 대립 유전자와 함께 4개 열을 모두 일치시켜 파일을 병합해야 합니다. 따라서 파일의 $1, $2, $4 및 $5 또는 $1, $2, $5 및 $4 열과 데이터 파일의 $1, $5, $6 및 $7 열이 정확하게 일치하는 항목을 찾아야 합니다. 데이터 파일에서 정확한 순서를 유지하는 것이 매우 중요합니다. 정렬할 수 없습니다(불행히도 조인을 사용하지 않습니다. 이는 이러한 유형의 질문에 대한 다른 인스턴스에서 찾은 제안 답변입니다).

나는 awk를 사용하여 수천 줄의 샘플 파일에서 작업할 수 있는 코드를 얻었지만 내 대규모 데이터 세트에 맞게 확장되지 않습니다(조회 파일에는 >3억 줄이 있고 데이터 파일에는 3천만 줄이 있음) - 아마도 코드 때문에 조회를 위해 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]
        }   
    }
}

관련 정보