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

debugcn 发表于 Dev

年龄87

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

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

（1）设置方法如下，尽管输出（groupVII 17978350）出现在列表中，但file_b的最后一行从输出中丢失了，尽管这会改变。我希望它按照设置的方式工作。

（2）如果名称具有特殊字符（例如冒号和连字符），则将其截断。我想在等号后加上全名。

（3）我想将file_b的条目/行与输出中的基因命中匹配，使得前两列是条目，第三列是基因命中。

file_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

file_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

output.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

尼科·甘布特（Niko Gambt）

# 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进行数据处理。

本文收集自互联网，转载请注明来源。

如有侵权，请联系[email protected] 删除。