AWK: частичное совпадение из одного файла в другом файле со специальным форматом в разных столбцах - PullRequest
0 голосов
/ 05 апреля 2020

Образец моего TSV file1 (который имеет больше дополнительных столбцов, но ради простоты уменьшен), где столбцы интересов CHROM и POS:

CHROM   POS         REF     ALT     QUAL    MoreColumns
chr11   8823729     G       C       605.77  ...
chr1    16619       C       T       95.77   ...
chr1    16949       A       C       559.77  ...
chr1    17005       A       G       172.77  ...
chr1    17020       G       A       345.77  ...
chr12   8822661     G       A       880.77  ...
chr1    17697       G       C       412.77  ...
chr14   8837474     T       C       411.77  ...
chr1    129285      G       A       2509.77 ...

A образец моего TSV file2, где интересующий столбец Extra_information и имеет следующий формат:

Column1 ... Column9     Extra_information                                                       Column11
data    ... longline    hg38:Chr12:8822661, hg19:Chr12:8975257, COM:morewords, dbSNP:link       No
data2   ... longline2   hg38:Chr11:8823729, hg19:chr12:8976325, COM:morewords2, dbSNP:link2     No
data3   ... longline3   hg38:chr12:8823762, hg19:Chr12:8976358, COM:morewords3                  Yes
data4   ... longline4   hg38:chr12:8835642, hg19:Chr12:8988238, dbSNP:link3                     No
data5   ... longline5   hg38:Chr14:8837474, hg19:chr12:8990070, dbSNP:link4                     Yes
data6   ... longline6   hg19:Chr12:8990937, COM:morewords4, dbSNP:link5                         No
data7   ... longline7   hg38:chr12:8839209, PC:someinfo                                         No

Моя проблема:

Я хочу выполнить частичное совпадение hg38:CHROM:POS от file1 до Extra_information от file2 и выведите строку file1 + "\t" 1, если частичное совпадение истинно, иначе строка file1 + "\t" 0. Chr также может быть chr в Extra_information из file2.

Мой желаемый первый выход

CHROM   POS         REF     ALT     QUAL    MoreColumns     Match
chr11   8823729     G       C       605.77  ...             1
chr1    16619       C       T       95.77   ...             0
chr1    16949       A       C       559.77  ...             0
chr1    17005       A       G       172.77  ...             0
chr1    17020       G       A       345.77  ...             0
chr12   8822661     G       A       880.77  ...             1
chr1    17697       G       C       412.77  ...             0
chr14   8837474     T       C       411.77  ...             1
chr1    129285      G       A       2509.77 ...             0

Мой предпочтительный второй выход

CHROM   POS         REF     ALT     QUAL    MoreColumns     Column1 ... Column9     Extra_information                                                       Column11
chr11   8823729     G       C       605.77  ...             data2   ... longline2   hg38:Chr11:8823729, hg19:chr12:8976325, COM:morewords2, dbSNP:link2     No
chr1    16619       C       T       95.77   ...             -       ... -               -                                                                   -
chr1    16949       A       C       559.77  ...             -       ... -               -                                                                   -
chr1    17005       A       G       172.77  ...             -       ... -               -                                                                   -
chr1    17020       G       A       345.77  ...             -       ... -               -                                                                   -
chr12   8822661     G       A       880.77  ...             data    ... longline    hg38:Chr12:8822661, hg19:Chr12:8975257, COM:morewords, dbSNP:link       No
chr1    17697       G       C       412.77  ...             -       ... -               -                                                                   -
chr14   8837474     T       C       411.77  ...             data5   ... longline5   hg38:Chr14:8837474, hg19:chr12:8990070, dbSNP:link4                     Yes
chr1    129285      G       A       2509.77 ...             -       ... -               -                                                                   -

Я попытался:

awk -F $'\t' 'NR == FNR 
    { a=("hg38:"file1[$1]":"file1[$2]); a=$NF; next } 
    { if ($10~$NF) {
       print file1[$0] "\t1"
    } else { 
       print file1[$0] "\t0"
    }
}' file1 file2

Как мне достичь желаемого результата (желательно второго), используя awk? (или если бы вы могли предложить любое другое решение bash)

Заранее спасибо.

Примечание: у меня ~ 70k строк из file1, чтобы выполнить частичное совпадение с file2 который содержит ~ 160 тыс. строк.

Редактировать:

По запросу @Hai Vu для полной строки:

File1:

https://drive.google.com/open?id=1kB4i7bpbA6zV1kRvGB3cBvt5RWYSurVJ

и File2:

https://drive.google.com/open?id=1gZ6qkYRuyEVT4Txom0sAawT2-F81reQN

1 Ответ

1 голос
/ 05 апреля 2020

Вот один из способов ее решения. Я создал скрипт AWK и назвал его hg38.awk . Чтобы вызвать его:

awk -f hg38.awk file2 file1

Обратите внимание, что я сканирую файл2 перед файлом1. Вот сценарий:

# In file2 where we found hg38
# We transform "hg38:Chr11:8823729," to "chr11:8823729"
# And use that as a key in the array `found`
NR == FNR && $4 ~ /^hg38:/ {
    extra = $4
    sub(/hg38:/, "", extra)
    sub(/Chr/, "chr", extra)
    sub(/,$/, "", extra)
    found[extra] = 1
}

# First line of file1
# Print the existing headers and an additional column
NR != FNR && FNR == 1 {
    print $0 "\tMatch"
    next
}

# Subsequent lines of file1
NR != FNR {
    printf $0
    key = $1 ":" $2
    if (key in found) {
        print "\t1"
    } else {
        print "\t0"
    }
}

Пояснения

  • В сценарии я сначала сканирую файл2 (см. Командную строку). Чтобы различать guish между двумя файлами, я смотрю на связь между NR и FNR переменными. Если они совпадают, я сканирую первый файл в командной строке (следовательно, file2). Если они отличаются, я работаю над файлом 1.
  • Для первого и второго блоков кода, я надеюсь, что комментарии будут адекватными
  • В последнем блоке я построил ключ из полей например, "chr11: 8823729" и проверьте, находится ли этот ключ в массиве found и выведите соответственно.

Update

Здесь я изменил скрипт для вывода вашего второго желаемый вывод. Изменения включают в себя сохранение всей строки файла2 и создание пустой строки.

NR == FNR && FNR == 1 {
    headers = $0
    empty_row = ""
    for (i = 0; i < NF; i++) {
        empty_row = "\t-" empty_row
    }
    next
}

# In file2 where we found hg38
# We transform "hg38:Chr11:8823729," to "chr11:8823729"
# And use that as a key in the array `found`
NR == FNR && $4 ~ /^hg38:/ {
    extra = $4
    sub(/hg38:/, "", extra)
    sub(/Chr/, "chr", extra)
    sub(/,$/, "", extra)
    found[extra] = $0
}

# First line of file1
# Print the existing headers and an additional column
NR != FNR && FNR == 1 {
    print $0 "\t" headers
    next
}

# Subsequent lines of file1
NR != FNR {
    printf $0
    key = $1 ":" $2
    if (key in found) {
        print "\t" found[key]
    } else {
        print empty_row
    }
}

Обновление 2

С последними данными от Google я обнаружил, что file2.tsv, поле 10 сложнее, чем я думал. После этого я смог разработать версию 3 своего решения:

# Works with TSV  (tab-separated values) file
BEGIN {
    FS = "\t"
}

# In file2.tsv, save the headers and create a row of empty data (just dashes)
NR == FNR && FNR == 1 {
    headers = $0
    empty_row = ""
    for (i = 0; i < NF; i++) {
        empty_row = "\t-" empty_row
    }
    next
}

# In file2.tsv where we found hg38
# We transform "hg38:Chr11:8823729," to "chr11:8823729"
# And use that as a key in the array `found`
NR == FNR && $10 ~ /^hg38:/ {
    extra = $10
    sub(/hg38:/, "", extra)
    sub(/Chr/, "chr", extra)
    sub(/,.*$/, "", extra)
    found[extra] = $0
}

# First line of file1
# Print the existing headers and additional columns
NR != FNR && FNR == 1 {
    print $0 "\t" headers
    next
}

# Subsequent lines of file1
NR != FNR {
    printf $0
    key = $1 ":" $2
    if (key in found) {
        print "\t" found[key]
    } else {
        print empty_row
    }
}
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...