题目来自生信技能树论坛
image.png
先从hg38的gtf中提取"ANXA1"基因
grep '"ANXA1"' hg38.gtf >ANXA1.gtf
在题目之前先分析要处理的数据的结构是什么样的以及要做什么
1.数据结构是怎么样的：
image.png
以gene开头，对应唯一的gene name
接下来是transcript 有唯一的transcript name
transcript下面的包括exon,CDS,start codon,stop codon,three_prime_utr,five_prime_utr都是属于transcript,都会对应有相同的transcript name
有的exon包括5'端的UTR区域或者3'端的UTR区域
而对于不包括UTR的exon就称之为CDS,所以CDS并不等同于exon
所以问题就很清楚了
2.要做什么：
如果要统计gene坐标：
则以gene name 为key建立gene的字典，方法是先构建dict_gene={}，然后直接把坐标存到dict_gene[gene_name]=[,]列表中储存坐标
如果要统计exon坐标：
则以gene name 为key建立exon的字典，然后字典中又以不同的transcript name 作为key建立字典，方法是先构建dict_exon[gene_name]={}，然后根据找到的transcript name,把transcript name作为key再存入，即dict_exon[gene_name][transcript_name]=[],把transcript对应的exon坐标存到列表中。
如果要统计UTR坐标：
同理；
image.png
image.png
python实现的代码如下
import collections
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib.patches as patches
matplotlib.style.use("ggplot")
gene_cord=collections.OrderedDict()
gene_exon=collections.OrderedDict()
gene_utr=collections.OrderedDict()
gene_start_codon=collections.OrderedDict()
gene_stop_codon=collections.OrderedDict()
with open ("ANXA1.gtf") as fh:
for line in fh:
lineL=line.strip().split("\t")
if lineL[2]=="gene":
gene_info=lineL[-1]
gene_infoL=gene_info.split("; ")
gene_name=gene_infoL[2]
gene_nameL=gene_name.split(" ")
gene_name=eval(gene_nameL[1])
gene_cord[gene_name]=[int(lineL[3]),int(lineL[4])] #对基因构建字典,key是基因名，value是基因坐标
gene_exon[gene_name]={}
gene_utr[gene_name]={}
elif lineL[2]=="transcript":
trans_info=lineL[-1]
trans_infoL=trans_info.split("; ")
trans_name=trans_infoL[9]
trans_nameL=trans_name.split(" ")
trans_name=eval(trans_nameL[1])
gene_exon[gene_name][trans_name]=[] #字典中套字典,最终的value先建立为空list
gene_exon[gene_name][trans_name].append(lineL[6]) #往空list中添加正负链信息的元素
gene_utr[gene_name][trans_name]=[]
elif lineL[2]=="exon":
gene_exon[gene_name][trans_name].extend([int(lineL[3]),int(lineL[4])]) #往list中加exon坐标
elif lineL[2] in ["five_prime_utr","three_prime_utr"]:
gene_utr[gene_name][trans_name].extend([int(lineL[3]),int(lineL[4])]) #往list中加UTR坐标
fig= plt.figure() #建一个画板
ax = fig.add_subplot(111) #将画板分为1行1列，图像画在从左到右从上到下第一块
trans_num=0
for trans_name,exon in gene_exon[gene_name].items():
trans_num+=1
if exon[0]=="+":
col="limegreen"
if exon[0]=="-":
col="magenta"
for i in range(1,len(exon),2): #注意！！！python中for i in range(n,m) 范围包括n,不包括m
#对每个exon画矩形图
#ax.add_patch(patches.Rectangle((x,y),width,height))
rect=patches.Rectangle((exon[i],trans_num-0.1),exon[i+1]-exon[i],0.2,color=col,fill=True)
ax.add_patch(rect)
#对每个intron画线图
if i < len(exon)-2:
#ax.plot([x1,x2],[y1,y2])
ax.plot([exon[i+1],exon[i+2]],[trans_num,trans_num],color="red",linewidth="0.5")
#对UTR画图
trans_num=0
for trans_name,utr in gene_utr[gene_name].items():
trans_num+=1
if utr==[]:
continue
else:
for i in range(0,len(utr),2):
rect=patches.Rectangle((utr[i],trans_num-0.1),utr[i+1]-utr[i],0.2,color="black")
ax.add_patch(rect)
ax.set_xlabel(gene_name,color="blue",fontsize=14,fontweight="bold")
ax.yaxis.set_major_locator(ticker.FixedLocator(range(1,trans_num+1)))
ax.set_yticklabels(gene_exon[gene_name].keys())
plt.xlim(gene_cord[gene_name])
plt.ylim([0.5,trans_num+0.5])
plt.show()
fig.savefig('rect.png', dpi=90, bbox_inches='tight')
几个注意的点：
1.list中增加元素的方法大概有两种方式：
第一种是list.append
另一种是list.expend
区别在于第一种是直接把整个模块加入list中，第二种是将元素一个个添加到list中
举个例子：
B = ['q', 'w', 'e', 'r']
B.append(['t', 'y'])
B
['q', 'w', 'e', 'r', ['t', 'y']]
A = ['q', 'w', 'e', 'r']
A.extend(['t', 'y'])
A
['q', 'w', 'e', 'r', 't', 'y']
我希望把exon和UTR坐标存到list中是当成每个元素去存的，方便后续处理 所以用的是extend
而trans_name是单个的，直接用append就好
2.for i in range(m,n,k),范围包括m，不包括n
网友评论