所以我有两个文件,一个 VCF 看起来像
88 Chr1 25 C - 3 2 1 1
88 Chr1 88 A T 7 2 1 1
88 Chr1 92 A C 16 4 1 1
还有一个基因看起来像
GENEID Start END
GENE_ID 11 155
GENE_ID 165 999
我想要一个脚本来查看第二个文件的第二个和第三个位置范围内是否存在基因位置(VCF 文件的第 3 列),然后将其打印出来。
到目前为止我所做的是加入文件并做
awk '{if (3>$12 && $3< $13) print }' > out
我所做的只是比较连接文件的当前行(它只在值在同一行时打印),我怎样才能让它比较第 3 列的所有行与第 12 和 13 列的所有行?
最好的, 塞尔格
最佳答案
我希望能有所帮助(编辑我更改代码以获得更高效的算法)
gawk '
#read input.genes and create list of limits (min, max)
NR == FNR {
#without header in input
if(NR>1) {
for(i=$2; i<=$3; i++){
limits[i]=limits[i]","$2"-"$3;
}
};
next
}
#read input.vcf, if column 3 is range of limits then print
{
if($3 in limits){
print $0, "between("limits[$3]")"
}
}' input.genes input.vcf
你得到:
88 Chr1 25 C - 3 2 1 1 between(,11-155)
88 Chr1 88 A T 7 2 1 1 between(,11-155)
88 Chr1 92 A C 16 4 1 1 between(,11-155)
python 中的这个算法针对使用字典的超大文件进行了优化
limits = [line.strip().split() for line in open("input.genes")]
limits.pop(0) #remove the header
limits = [map(int,v[1:]) for v in limits]
dict_limits = {}
for start, finish in limits:
for i in xrange(start, finish+1):
if i not in dict_limits:
dict_limits[i] = []
dict_limits[i].append((start,finish))
OUTPUT = open("my_output.txt", "w")
for reg in open("input.vcf"):
v_reg = reg.strip().split()
if int(v_reg[2]) in dict_limits:
OUTPUT.write(reg.strip() + "\tbetween({})\n".format(str(dict_limits[int(v_reg[2])])))
OUTPUT.close()
你得到:
88 Chr1 25 C - 3 2 1 1 between([(11, 155)]) 88 Chr1 88 A T 7 2 1 1 between([(11, 155)]) 88 Chr1 92 A C 16 4 1 1 between([(11, 155)])
关于python - 将一列值与linux环境中的所有列进行比较,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/30191109/