搜索
楼主: Jimmy

生信编程直播第五题:根据GTF画基因的多个转录本结构

[复制链接]

0

主题

4

帖子

41

积分

新手上路

Rank: 1

积分
41
发表于 2017-3-12 23:30:35 | 显示全部楼层
这道题确实太难了,花了2个晚上4个小时的时间,才做出来perl+R,图片还需要一些美观。不过差不多了!
perl部分
#!/usr/bin/perl
use strict;
my ($gene_name,$row,$chr,$gene_start,$gene_end,@cells,$x,$count);
open INPUT,"Btau_5.0.1_genomic.gff";
open OUTPUT,">graph_data.txt";
my $trans_num=0;
print OUTPUT "transcript_ID\ttrans_num\tpos\n";
while(<INPUT>){
        chomp;
        if ((/transcript_id=XM/)&&(/gene=ABI3BP;/)){
                @cells=split/\t/,$_;
                if($cells[2] eq "mRNA"){
                        $trans_num++;
                }
                if($cells[2] eq "exon"){
                        /transcript_id=(XM.*?);/;
                        $gene_start=$cells[3];
                        $gene_end=$cells[4];
                        my $transcript_ID=$1;
                        for ($x=$cells[3];$x<=$cells[4];$x += 1){
                                print OUTPUT "$transcript_ID\t$trans_num\t$x\n";
                        }
                }
        }
}
R部分
library(ggplot2)
a<-read.table("D:/R_figure/script_practice/graph_data.txt",header=T)
p<-ggplot(a,aes(pos,trans_num))+geom_point(col="green")+theme_set(theme_bw())+theme(panel.grid.minor=element_blank())+theme(panel.grid.major=element_blank())+theme(panel.background=element_blank())
p<-p+scale_y_continuous(breaks=seq(1,32,1),labels=c("XM_015461955.1","XM_010824537.2","XM_010824570.2","XM_015461982.1","XM_015461965.1","XM_015461963.1","XM_015461981.1","XM_015461978.1","XM_015461969.1","XM_015461967.1","XM_015461964.1","XM_015461961.1","XM_015461962.1","XM_015461980.1","XM_015461979.1","XM_015461977.1","XM_015461976.1","XM_015461956.1","XM_015461975.1","XM_015461966.1","XM_015461974.1","XM_015461968.1","XM_015461959.1","XM_015461973.1","XM_015461972.1","XM_015461971.1","XM_015461958.1","XM_015461970.1","XM_015461960.1","XM_015461957.1","XM_015461954.1","XM_010824535.2"))+labs(x="chr1(bp)",y="transcript ID")
回复 支持 反对

使用道具 举报

0

主题

19

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-3-13 05:06:33 | 显示全部楼层
本帖最后由 朝曦曦 于 2017-3-13 19:16 编辑

学员 205 Python and perl

总算刚写出来第一部分,我用的是Homo_sapiens.GRCh38.87.chr.gtf 作为模板,然而没有找到UTR 区域在哪里。别的和讲师给的答案是一致的。
[AppleScript] 纯文本查看 复制代码
#### Parse exon coordinates of each transcript
from collections import OrderedDict

geneCord = OrderedDict()
GRCh38Exon = OrderedDict()
GRCh38UTR = OrderedDict()

with open('Homo_sapiens.GRCh38.87.chr.gtf', 'rt') as f:
    for line in f:
        line = line.rstrip()
        if line.startswith('#'):
            continue
        
        lst = line.split('\t')
        
        if lst[2].startswith('gene'):
            geneInfo = lst[-1].split(';')
            geneName=geneInfo[2].split('"')[1]
            geneCord[geneName] = [int(lst[3]),int(lst[4])]
            GRCh38Exon[geneName] = OrderedDict()
            GRCh38UTR[geneName] = OrderedDict()
        
        elif lst[2].startswith('transcript'):
            trptInfo = lst[-1].split(';')
            trptName = trptInfo[9].split('"')[1]
            GRCh38Exon[geneName][trptName] = []
            GRCh38Exon[geneName][trptName].append(lst[6])
            GRCh38UTR[geneName][trptName] = []
            
        elif lst[2] == 'exon':
            GRCh38Exon[geneName][trptName].extend([int(lst[3]), int(lst[4])])
答案长这样,
[AppleScript] 纯文本查看 复制代码
OrderedDict([('ANXA1-002',
              ['+',
               73151757,
               73151924,
               73158522,
               73158601,
               73158695,
               73158803,
               73159329,
               73159423,
               73160263,
               73160376,
               73160803,
               73160893,
               73162782,
               73162861,
               73163476,
               73163532,
图也是没有UTR区域的,如下:

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 反对

使用道具 举报

0

主题

6

帖子

67

积分

注册会员

Rank: 2

积分
67
发表于 2017-3-29 11:18:45 | 显示全部楼层
本帖最后由 highwind 于 2017-3-29 11:39 编辑

A:问题拆解1(看看数据)
之前大神都说了gtf是该死的格式还是genepred好,可是转换格式的工具只能在unix下运行,暂时还是老老实实做题吧
看起来又是个数组的形式

问题拆解2(想想算法)
提取信息条目,然后用matplotlib作图
作图的时候可以用图层的方式一个个组件叠加上去

问题拆解3(算法实现)
[Python] 纯文本查看 复制代码
# -*- coding: utf-8 -*-
import time
starttime = time.clock()

def find_target_data(gene_name, gtf_file_path, chunk=2048*2048): 
    with open(gtf_file_path) as f:
        target_data = []
        get = 0 
        buffer_list = ['']
        print('Wait ...\n')
        while True:
            buffer = buffer_list[-1] + f.read(chunk)
            print(buffer[buffer.find('\n')+1:buffer.find('\n')+3], end='\r')
            if ('gene_name "'+gene_name.upper()+'"') in buffer:
                buffer_list = buffer.split('\n')
                for line in buffer_list[:-1]:
                    if ('gene_name "'+gene_name.upper()+'"') in line:
                        get = 1
                        line_list = line.split('\t')
                        target_data.append(line_list)

            else:
                if get == 1:
                    break
            if len(buffer) < chunk:
                break #这里不会把小于chunk的部分给忽略掉吗?
        if not target_data:
            print('wrong name\n')
            raise NameError('Gene name not correct\n')#这个学习了,牛逼闪闪的错误提示
        print('Gene Get\n')
    return(target_data)#这是一个二维数组
        
def plot_gene(target_data, gene_name, plot_path='', line_width=5):
    import matplotlib.pyplot as plt
    import matplotlib.ticker as ticker
    import re
    i = 0
    j = 0
    transcript_list = ['']
    Fig = plt.figure()
    ax = Fig.add_subplot(111)#111说是1行1列的第1幅图的意思
    color_dict = {'gene':'black',\
                  'transcript':'blue',\
                  'CDS':'yellow',\
                  'exon':'cyan',\
                  'start_codon':'red',\
                  'stop_codon':'magenta',\
                  'UTR':'green'}
    while(i < len(target_data)):#用图层的方法叠加,但是叠加的层次倒还需要算法调整,下面就偷懒了,只用了文件里默认的顺序
        START = int(target_data[i][3])
        END = int(target_data[i][4])
        if(target_data[i][2] == 'gene'):
            g_START = START
            g_END = END
        elif(target_data[i][2] == 'transcript'):
            transcript_name = re.findall('transcript_name "(.*?)"',target_data[i][8])[0]
            transcript_list.append(transcript_name)
            j += 1
            ax.plot([g_START,g_END],[j,j],linewidth=1,color=color_dict['gene']) 
            if(target_data[i][6] == '+'):
                ax.plot(g_END,j,">",color=color_dict['gene']) #没想到画个箭头居然还挺容易
            else: #应该相当于elif(target_data[i][6] == '-')
                ax.plot(g_START,j,"<",color=color_dict['gene'])
            ax.plot([START,END],[j,j],linewidth=4.0,color=color_dict['transcript']) 
        elif(target_data[i][2] == 'UTR'):
            ax.plot([START,END],[j,j],linewidth=8.0,color=color_dict['UTR'])
        elif(target_data[i][2] == 'CDS'):
            ax.plot([START,END],[j,j],linewidth=8.0,color=color_dict['CDS'])
        elif(target_data[i][2] == 'exon'):
            ax.plot([START,END],[j,j],linewidth=8.0,color=color_dict['exon'])
        elif(target_data[i][2] == 'start_codon'):
            ax.plot([START,END],[j,j],'d',color=color_dict['start_codon']) #d表示菱形
        elif(target_data[i][2] == 'stop_codon'):
            ax.plot([START,END],[j,j],'d',color=color_dict['stop_codon']) #d表示菱形
        else:
            pass
        i += 1

    #调整x轴标签用的一些计算
    import numpy as np #学习一下使用np中的间隔数函数linspace
    t_START = int(g_START/100-1)*100
    t_END = int(g_END/100+1)*100
    n=5
    ax.xaxis.set_major_locator(ticker.FixedLocator(np.linspace(t_START,t_END,n))) #这里用ax.set_xticklabels(range(t_START,t_END))显示是不对的
    ax.ticklabel_format(style='plain') #以为应该用ax.ticklabel_format(useOffset=False)的,但是发现这个似乎是改y轴科学记数法的?
    ax.tick_params(left='off')#干掉左边tick
    ax.set_ylabel('Trascripts',fontsize=16)
    ax.set_yticklabels(transcript_list)
    ax.yaxis.set_major_locator(ticker.FixedLocator(np.linspace(0,len(transcript_list),len(transcript_list)+1)))#为了调整转录本少的基因可真是不容易
    ax.spines['top'].set_visible(False)   #上边框不可见
    ax.spines['right'].set_visible(False) #右边框不可见
    ax.spines['left'].set_visible(False)  #左边框不可见
    ax.set_title('Gene Structure of '+gene_name.upper(),fontsize=16)
    #画图例
    import matplotlib.patches as mpatches
    patch_list = []
    for c in color_dict:
        patch_list.append(mpatches.Patch(color=color_dict[c],label=c))
    ax.legend(handles=patch_list,\
              bbox_to_anchor=(1.35, 0),loc=4)  #这里的1.35是完全手动调节出来的,不知道大家还有什么别的经验?  
    Fig.savefig(gene_name,bbox_inches='tight', pad_inches=1) #为了使得fig.show()和savefig结果一致还是要做些功课的

gtf_file_path = 'gencode.v7.annotation_goodContig.gtf'
gene_name = 'egfr'
target_data = find_target_data(gene_name, gtf_file_path)
plot_gene(target_data, gene_name)
endtime=time.clock()
print("本次计算用时{:.3f}秒".format(endtime - starttime)) 

计算结果:
本次计算用时6.522秒

不过有些只有比较少转录本的基因的图看起来就比较丑了:


而且这个结果里有个误解:那个标识为蓝色的“transcript”并不是真的transcript,因为这里显示的长度是在基因组上而非转录组上的!
也就是说这段蓝色只是说明了在基因组上转录本的起始范围,而不是真的transcript就这么长的一段,
要不然一个start codon空了好一段才接上CDS也说不过去啊。

问题拆解4(算法优化)
这个。。我真的没太看懂Buffer那部分的buffer = buffer_list[-1] + f.read(chunk)。还有这个chunk怎么优化呢?

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 1 反对 0

使用道具 举报

0

主题

11

帖子

71

积分

注册会员

Rank: 2

积分
71
发表于 2017-4-8 18:43:36 | 显示全部楼层
x2yline 发表于 2017-2-15 19:37
77-R-python-shell x2yline

数据来源:(第三题的文件)ftp://ftp.ensembl.org/pub/release-87/gtf/homo_s ...

这个代码写的太经典了,涉及的chunk,buffer相关的内容,都是以前未接触的,不知道哪里有好的解释。
回复 支持 反对

使用道具 举报

2

主题

34

帖子

773

积分

高级会员

Rank: 4

积分
773
发表于 2017-4-8 23:13:21 | 显示全部楼层
本帖最后由 x2yline 于 2017-4-8 23:22 编辑
13235695523 发表于 2017-4-8 18:43
这个代码写的太经典了,涉及的chunk,buffer相关的内容,都是以前未接触的,不知道哪里有好的解释。 ...

看dongye老师第二题的python视频就知道buffer这个东西了,很简单的chunk值就是字节大小,chunk如果为1那f.read(chunk)则只读取一个字节,每次的buffer=f.read(chunk)就只读取了一个字母, 如果chunk=2048*2048那一次就读取2048*2048=4194304个字节,每次的buffer=f.read(chunk)就读取了4194304个字母.

我也是看了dongye老师第二题的视频, 才知道可以这样来加速

回复 支持 反对

使用道具 举报

2

主题

34

帖子

773

积分

高级会员

Rank: 4

积分
773
发表于 2017-4-8 23:51:20 | 显示全部楼层
本帖最后由 x2yline 于 2017-4-9 07:42 编辑
highwind 发表于 2017-3-29 11:18
A:问题拆解1(看看数据)
之前大神都说了gtf是该死的格式还是genepred好,可是转换格式的工具只能在unix下 ...

建议看看第二题视频,阅读一下f.read()函数的使用说明, 多试一试f.read()这个函数,基本就懂buffer了
回复 支持 反对

使用道具 举报

1

主题

4

帖子

52

积分

注册会员

Rank: 2

积分
52
发表于 2017-8-17 14:32:15 | 显示全部楼层
这个帖子超赞,手动点赞搜藏,找个时间再来撸一遍
回复 支持 反对

使用道具 举报

4

主题

48

帖子

778

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
778
发表于 2017-10-3 09:11:08 | 显示全部楼层
x2yline 发表于 2017-2-15 19:37
77-R-python-shell x2yline

数据来源:(第三题的文件)ftp://ftp.ensembl.org/pub/release-87/gtf/homo_s ...

你可以试试data.table包  这个c写的,读取1G以上文本,速度比read.table等快太多。另外一个readr包也可以读取大数据,比data.table包慢一点。data.table有个缺陷,没有row.names参数。这个读取,进入需要自己调整打到自己满意的格式。另外一个注意的就是data.table的fread读取一个文本,默认是data.table 和 data.frame,在读取的时候里面加一个参数。
回复 支持 反对

使用道具 举报

1

主题

7

帖子

164

积分

注册会员

Rank: 2

积分
164
发表于 2017-10-4 00:41:45 | 显示全部楼层
本帖最后由 源氏 于 2017-10-4 11:37 编辑

我来了!!
这道题目我终于有机会来回答了
25行代码, 从提取到画图, 不用任何转换!!直接提取

[Python] 纯文本查看 复制代码
!grep '"ANXA1"' Homo_sapiens.GRCh38.90.chr.gtf|awk '{print$3"\t"$4-1"\t"$5-1}'>data2
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
plt.xkcd()
track = 0
plt.figure(figsize=(12,6))
color_dic={'exon':'#00896C','CDS':'#78552B', 'start_codon':'#CB1B45','stop_codon':'black','five_prime_utr':'#F19483','three_prime_utr':'#2EA9DF'}
linewith_dic={'exon':8.0,'CDS':6.0, 'start_codon':8.0,'stop_codon':8.0,'five_prime_utr':4.0,'three_prime_utr':4.0}
with open('data2','r') as fr:
    lines = fr.readlines()
    gene_start = lines[0].rstrip().split()[1]
    gene_end = lines[0].rstrip().split()[2] 
    plt.xlim(int(gene_start) - 500,int(gene_end)+ 500)
    del lines[0]
    for line in lines:
        line = line.rstrip().split()
        if line[0] == 'transcript':
            track += 1
            plt.plot([int(line[1]),int(line[2])],[track,track],color="black")
        else:
            plt.plot([int(line[1]),int(line[2])],[track,track],color=color_dic[line[0]],linewidth=linewith_dic[line[0]])
    plt.title("the transcripts of  XXX  gene")
    plt.yticks(np.arange(0,track+3))
plt.show();


效果图:

只要修改基因名,就可以完成从GTF到最后结果图
其实GTF数据很好提取
TP53的如图:

夜深, 图例什么就先不弄了,
晚安,生信人!!

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复 支持 2 反对 0

使用道具 举报

5

主题

32

帖子

484

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
484
发表于 2017-10-7 16:42:52 | 显示全部楼层
突然意识到我们这里指的基因的范围并没有包括启动子,增强子序列,因为他们不会被转录,他们是在转录起始位点前面
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-21 01:14 , Processed in 0.036998 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.