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 collectionsimport matplotlibimport matplotlib.pyplot as pltimport matplotlib.ticker as tickerimport 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=0for 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

ax.plot([exon[i+1],exon[i+2]],[trans_num,trans_num],color="red",linewidth="0.5")#对UTR画图            trans_num=0for 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

作者：天秤座的机器狗

链接：https://www.jianshu.com/p/abb29ce1f76d