主页 / unix / 问题 / 495683
Accepted
Age87
Asked: 2019-01-21 17:08:29 +0800 CST

使用 file_b 中 2 列的信息从 file_a 中提取名称

  • 772

建立在：从 File_B 中提取与 File_A 有重叠间隔的名称

我想通过将file_b的第一列与file_a匹配来提取基因名称（通常是第10列，“Name =”之后的内容），如果file_b的第二列位于基因区间内，则提取基因名称，由第4列和file_b 的 5 个。第一列必须匹配，这样我每行只能得到一个基因（file_b），但理论上我可以让多个相邻行（column_b）匹配同一个基因（例如，如果 file_b 中的第二行是MT 4065）

我现在的代码有一些问题。

(1) 按照下面的设置方式，file_b 的最后一行从输出中丢失，尽管如果此行 ( groupVII 17978350) 出现在列表中，情况会发生变化。我希望它按照设置的方式工作。

(2) 如果名称包含特殊字符（例如冒号和连字符），则名称将被截断。我想在等号后面加上全名。

(3) 我想将 file_b 的条目/行与输出中的基因匹配匹配，这样前两列是条目，第三列是基因匹配。

文件_a.tsv

MT  insdc   gene    2851    3825    .   +   .   ID=gene:ENSGACG00000020925  Name=mt-nd1 biotype=protein_coding  description=NADH dehydrogenase 1%2C mitochondrial [Source:ZFIN%3BAcc:ZDB-GENE-011205-7] gene_id=ENSGACG00000020925  logic_name=mt_genbank_import    version=1
MT  insdc   gene    4036    5082    .   +   .   ID=gene:ENSGACG00000020929  Name=mt-nd2 biotype=protein_coding  description=NADH dehydrogenase 2%2C mitochondrial [Source:ZFIN%3BAcc:ZDB-GENE-011205-8] gene_id=ENSGACG00000020929  logic_name=mt_genbank_import    version=1
groupIII    ensembl gene    7332324 7334769 .   -   .   ID=gene:ENSGACG00000015265  Name=si:dkeyp-68b7.10   biotype=protein_coding  description=si:dkeyp-68b7.10 [Source:ZFIN%3BAcc:ZDB-GENE-070912-667]    gene_id=ENSGACG00000015265  logic_name=ensembl  version=1
groupIV ensembl gene    1368026 1374881 .   +   .   ID=gene:ENSGACG00000016447  Name=hnrnpa0b   biotype=protein_coding  description=heterogeneous nuclear ribonucleoprotein A0b [Source:ZFIN%3BAcc:ZDB-GENE-030131-6154]    gene_id=ENSGACG00000016447  logic_name=ensembl  version=1
groupIV ensembl gene    5347339 5349041 .   -   .   ID=gene:ENSGACG00000017010  Name=zgc:153018 biotype=protein_coding  description=zgc:153018 [Source:ZFIN%3BAcc:ZDB-GENE-060929-752]  gene_id=ENSGACG00000017010  logic_name=ensembl  version=1
groupV  ensembl gene    120615  125489  .   +   .   ID=gene:ENSGACG00000002103  Name=zdhhc6 biotype=protein_coding  description=zinc finger%2C DHHC-type containing 6 [Source:ZFIN%3BAcc:ZDB-GENE-030131-3189]  gene_id=ENSGACG00000002103  logic_name=ensembl  version=1
groupVI ensembl gene    11230354    11232784    .   +   .   ID=gene:ENSGACG00000009527  Name=bnip4  biotype=protein_coding  description=BCL2 interacting protein 4 [Source:ZFIN%3BAcc:ZDB-GENE-051113-212]  gene_id=ENSGACG00000009527  logic_name=ensembl  version=1
groupVII    ensembl gene    2271611 2277214 .   +   .   ID=gene:ENSGACG00000019012  Name=sf3b2  biotype=protein_coding  description=splicing factor 3b%2C subunit 2 [Source:ZFIN%3BAcc:ZDB-GENE-070928-1]   gene_id=ENSGACG00000019012  logic_name=ensembl  version=2
groupVII    ensembl gene    15815857    15824549    .   +   .   ID=gene:ENSGACG00000020296  Name=mpp1   biotype=protein_coding  description=membrane protein%2C palmitoylated 1 [Source:ZFIN%3BAcc:ZDB-GENE-031113-4]   gene_id=ENSGACG00000020296  logic_name=ensembl  version=1
groupVII    ensembl gene    17978322    17982388    .   +   .   ID=gene:ENSGACG00000020399  Name=si:ch211-284e13.4  biotype=protein_coding  description=si:ch211-284e13.4 [Source:ZFIN%3BAcc:ZDB-GENE-060526-161]   gene_id=ENSGACG00000020399  logic_name=ensembl  version=1

文件_b.tsv

MT  4050
groupIII    7332350
groupIV 5347350
groupVI 11230375
groupVII    17978350

代码：

while read -r id pos; do awk -v id="$id" -v pos="$pos" '$1 == id && pos > $4 && pos < $5 { if (gensub(/.*Name=([A-Za-z0-9]*).*/, "\\1", 1) !~ /\s/) print gensub(/.*Name=([A-Za-z0-9]*).*/, "\\1", 1); }' <file_a.tsv; done < file_b.tsv > output.tsv

输出.tsv

mt
si
zgc
bnip4

所需输出

MT  4050    mt-nd2
groupIII    7332350 si:dkeyp-68b7.10
groupIV 5347350 zgc:153018
groupVI 11230375    bnip4
groupVII    17978350    si:ch211-284e13.4
bash shell-script

2 个回答

  1. 您预期的显示输出在我看来不一致（2 行--> 第一和第三行），如果这是一个错字，那么您能否尝试以下操作。

    awk 'FNR==NR{a[$1]=$2;next} ($1 in a) && (a[$1]>=$4 && a[$1]<=$5){sub("Name=","",$10);print $1,a[$1],$10}'  b.tsv a.tsv > output.tsv
    • 1
  2. Best Answer
    # save this as script.awk or whatevernameyouwant.awk

function within_range(val, lower, upper, proximity) {
    # you can specify the "proximity" as required
    return val > lower - proximity && val < upper + proximity
}

BEGIN {
    OFS="\t"
}

$1 == id && within_range(pos, $4, $5, 100) {
    name = gensub(/.*Name=([^\t]*).*/, "\\1", 1)
    if (name ~ /[^[:space:]]+/)
        print id, pos, name
}

    然后运行

    while read -r id pos
do
    awk -v id=$id -v pos=$pos -f script.awk file_a.tsv
done < file_b.tsv > output.tsv

    请确保.tsv文件中的字段在处理它们之前用制表符分隔。我的输出：

    MT  4050    mt-nd2
groupIII    7332350 si:dkeyp-68b7.10
groupIV 5347350 zgc:153018
groupVI 11230375    bnip4
groupVII    17978350    si:ch211-284e13.4

    对于ID MT，基因命中应该mt-nd2没有mt-nd1。

    我仍然推荐使用 Python 进行数据处理。

    • 1

相关问题