生信编程实战第5题（python）

题目来自生信技能树论坛

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