|
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享： 问题需求的描述- 任意给定基因组的chr：pos，判断它在那个基因上面？
基因的chr，start，end都是已知的（可以从UCSC或者NCBI获得）- 已知CNV数据在指定染色体上的chr：2075000-2930999，怎样批量获取其对应的基因？
- 如果该位置位于基因上，那么是位于内含子还是外显子，进一步说是第几个外显子或者第几个内含子
- 如果不在基因上，那么指定坐标距离哪一个基因的TSS（转录起始位点）最近？
以上这些问题总结起来就是对指定位点或者空间进行genomic features注释，针对这个问题呢，已经有annovar，snpEFF，VEP可以进行基因组特征的注释。 简单 的练习为判断指定参考基因组坐标位置是否在基因上，在哪个基因上？ 以下的python版本是参照曾老师的perl版本（http://www./thread-1321-1-1.html）写的 #!/usr/bin/python
from collections import OrderedDict
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('my_pos',help='input the position',type=int)
args = parser.parse_args()
my_pos=args.my_pos
#定义一个函数找出最邻近的位点
def nearest_search(array,value,low=0,high=0):
if high==0:
high = len(array)-1
else:
high = high
if low==0:
low = 0
else:
low = low
#print(low,high)
if array[low] > array[high]:
print("Array isn't sorted")
else:
if value < array[low] :
#print(low,array[low])
if value < (array[low]+array[low-1])/2:
return(low-1)
else:
return(low)
elif value > array[high]:
#print(high,array[high])
if value < (array[high]+array[high-1])/2:
#print(high)
return(high)
else:
#print(high+1)
return(high+1)
mid = int((low + high)/2)
if array[mid] == value:
return(mid)
elif array[mid] < value:
return nearest_search(array,value,mid+1,high)
else:
return nearest_search(array,value,low,mid-1)
#创建一个字典，以染色体号为键，每个基因的转录起始位点为值
Chr=OrderedDict()
#创建一个字典，以染色体号和转录起始位位点为键，基因的信息为值
gene_info=OrderedDict()
with open('hg38.tss','r') as f:
for line in f:
arr = line.strip().split('\t')
key=arr[1]
if key not in Chr:
Chr[key]=[]
Chr[key].append(int(arr[2]))
else:
Chr[key].append(int(arr[2]))
key2=arr[1]+'\t'+arr[2]
gene_info[key2]=line.strip()
my_list = sorted(Chr['chr1'])
nearest_index=nearest_search(my_list,my_pos)
print(my_list[nearest_index])
key='chr1'+'\t'+str(my_list[nearest_index])
print(key)
print(gene_info[key])
arr = gene_info[key].split('\t')
gene=arr[0]
if my_pos < my_list[nearest_index]:
print('%s is nearest gene of %s' % (gene,my_pos))
else:
print('%s within %s' % (my_pos,gene))
结果文件： $ python anno_1.py 27371
27370
chr1 27370
NR_024540 chr1 27370 31370 1
27371 within NR_024540
$ python anno_1.py 27366
27370
chr1 27370
NR_024540 chr1 27370 31370 1
NR_024540 is nearest gene of 27366
文末友情推荐
|
