所以我用awk/bash写代码。我有两个文件，第一个文件的格式如下：

chr1_2376428_A_T
chr1_5465765_T_A
chr1_8958392_C_G
....
chrM_237426_C_G

字符串
这个文件跨越所有染色体，每个染色体有多个样本。现在我的第二个文件是一个基因注解文件，其格式如下：

chr1 HAVANA Exon 13642 13859 ....
chr1 HAVANA START_CODON 13642 13859 ....
....

型
如果染色体匹配，我想搜索文件1中的值是否落在注解文件中开始（$4）和结束（$5）坐标所描述的任何区域之间。
示例文本文件：

echo -e chr1_2376428_A_T"/n"chr1_5465765_T_A"/n"chr1_8958392_C_G"/n"chr1_9542312_T_G > coordinates.txt

echo -e chr1"\t"HAVANA"\t"Exon"\t"13642"\t"13859"\n"chr1"\t"HAVANA"\t"START_CODON"\t"13642"\t"13859"\n"chr1"\t"HAVANA"\t"Exon"\t"5465750"\t"5465810"\n"chr1"\t"HAVANA"\t"gene"\t"6465750"\t"6465810"\n"chr1"\t"HAVANA"\t"Exon"\t"8958350"\t"8958461"\n" | gzip -c > gencode.v34.annotation.gtf.gz

型
我使用了以下代码：

zcat gencode.v34.annotation.gtf.gz | awk -v File=/pathto/Coordinates.txt 'BEGIN{comm="cat "File; while(comm|getline){split($1,a,"_");dict[a[1]]=a[2];}}{if ($1 in dict){if (dict[$1]>$4 && dict[1]<$5){print $1"_"dict[$1]"_"$3}}}'

型
这里的问题是，我得到的每个染色体的结果是每个染色体的最后一个可能的样本，这使我认为只有键的最后一个值（染色体编号）被比较和/或只被记录。
目前的成果：

None

型
预期成果：

chr1_5465765_Exon
chr1_8958392_Exon

型
谢谢你
编辑：下面是上面的awk代码，gawk -o-打印得很好，很清晰：

BEGIN {
        comm = "cat " File
        while (comm | getline) {
                split($1, a, "_")
                dict[a[1]] = a[2]
        }
}

{
        if ($1 in dict) {
                if (dict[$1] > $4 && dict[1] < $5) {
                        print $1 "_" dict[$1] "_" $3
                }
        }
}

型

$ cat >sqlite-commands <<'EOD' create table a (chr, v2, tag, start, end); .separator \t .import annotations_file a create table b (chr, idx, v3, v4); .separator _ .import samples_file b select b.chr, b.idx, a.tag from a join b on a.chr = b.chr where a.start < b.idx and b.idx < a.end; EOD $ cat -A annotations_file chr1^IHAVANA^IExon^I13642^I13859$ chr1^IHAVANA^ISTART_CODON^I13642^I13859$ chr1^IHAVANA^IExon^I5465750^I5465810$ chr1^IHAVANA^Igene^I6465750^I6465810$ chr1^IHAVANA^IExon^I8958350^I8958461$ $ cat samples_file chr1_2376428_A_T chr1_5465765_T_A chr1_8958392_C_G chr1_9542312_T_G $

1条答案

按热度按时间

ohtdti5x1#

array[key]=value只允许每个键有一个值-每当设置一个新值时，它就会被覆盖（即每个dict[a[1]] = a[2]）
dict[1]似乎是脚本后面的一个错字

使用生物信息学特定的工具，如bedtools，正如@timur-shtatland在评论中提到的，是首选。然而，使用SQLite的朴素SQL解决方案可能是：

字符串
然后又道：

$ sqlite3 '' <sqlite-commands
chr1_5465765_Exon
chr1_8958392_Exon
$

型
我理解传递''是因为数据库使用临时文件，并允许处理内存中无法容纳的数据。
要从压缩文件加载数据而不解压到磁盘，可以使用fifo：

$ mkfifo annotations_fifo
$ zcat annotations_file.gz >annotations_fifo &

型
并在sqlite import元命令中使用fifo的名称（即annotations_fifo）。（记住之后删除fifo。）
或者如果使用bash和Linux，也许可以使用进程重定向：

$ sqlite3 '' <sqlite-commands 4< <(zcat annotations_file.gz)

型
并在import元命令中使用/proc/self/fd/4作为文件名：

赞(0）回复(0）举报 12个月前

unix 如何在awk中访问字典中一个键的多个值？

1条答案

相关问题

热门标签

最新问答