搜索
楼主: Jimmy

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

[复制链接]

2

主题

34

帖子

777

积分

高级会员

Rank: 4

积分
777
发表于 2017-2-16 00:12:04 | 显示全部楼层
Jimmy 发表于 2017-2-1 09:10
聪明,这个代码是我博客历史遗留问题,你去看我博客就明白了!

legend其实可以放在作图区域外侧的http://stackoverflow.com/questio ... ea-in-base-graphics
回复 支持 反对

使用道具 举报

0

主题

15

帖子

151

积分

注册会员

Rank: 2

积分
151
发表于 2017-2-16 11:13:43 | 显示全部楼层
用perl提取出来应该只是第一步啊 ?昨晚第一步就进行不下去了啊,意思是后面画图还是需要用python或者R嘛 ?
回复 支持 反对

使用道具 举报

0

主题

12

帖子

165

积分

注册会员

Rank: 2

积分
165
发表于 2017-2-17 10:33:13 | 显示全部楼层
107 Dorom
先是看了perl的视频,发现好像perl并不适合作图,而且视频貌似还没有录完的赶脚。然后看了R的视频,本来是菜鸟一只,自己根本就写不出,所以只有跟着视频里的一步一步完成,最后图是画好了,而且自己还换了个基因试了一下。每行命令还没有理解透,还得继续加油,多多练习
回复 支持 反对

使用道具 举报

6

主题

23

帖子

201

积分

中级会员

Rank: 3Rank: 3

积分
201
发表于 2017-2-17 23:23:54 | 显示全部楼层
本帖最后由 qin_qinyang 于 2017-2-17 23:25 编辑

R-123-qin
还是习惯先把数据拆分成行,然后一行一行运行,看一下每个命令运行的结果到底是什么,才理解的透彻一些。不过通过拆分脚本,发现数据的整体结构,好像也不是那么难。
[AppleScript] 纯文本查看 复制代码
#!/bin/python
from collections import OrderedDict
import re
geneCord = OrderedDict()
GRCh38Exon = OrderedDict()
GRCh38UTR = OrderedDict()

with open('/Users/yangqin/Desktop/biotree/python/lession_5/temp1.gtf', 'rt') as f:
        for line in f:
                line = line.rstrip()
                if line.startswith('#'):
                        continue
                lst = line.split('\t')
                if lst[1] =='havana':
                        if lst[2] == 'gene':
                                geneName = re.search(r'gene_name \"([^;]+)\";?',lst[-1]).group(1)
                                geneCord[geneName] = [int(lst[3]),int(lst[4])]
                                GRCh38Exon[geneName] = OrderedDict()
                                GRCh38UTR[geneName] = OrderedDict()
                        elif lst[2] == 'transcript':
                                trpName = re.search(r'transcript_id \"([^;]+)\";?', lst[-1]).group(1)
                                GRCh38Exon[geneName][trpName] = []
                                GRCh38Exon[geneName][trpName].append(lst[6])
                                GRCh38UTR[geneName][trpName] = []
                        elif lst[2] == 'exon':
                                GRCh38Exon[geneName][trpName].extend([int(lst[3]),int(lst[4])])
                        elif lst[2] in ['five_prime_utr','three_prime_utr']:
                                GRCh38UTR[geneName][trpName].extend([int(lst[3]),int(lst[4])])
# geneCord['DDX11L1']
# GRCh38Exon['DDX11L1']
# GRCh38UTR['DDX11L1']
for k,v in GRCh38Exon['DDX11L1'].items():
         print (k,v)
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib.patches as patches
from matplotlib.patches import Rectangle
matplotlib.style.use('ggplot')

fig = plt.figure()
ax = fig.add_subplot(111)

trpNum = 0
for trpt,exon in GRCh38Exon['DDX11L1'].items():
        trpNum += 1
        if exon[0] == '+':
                col ='blue'
        elif exon[0] == '-':
                col = 'red'
        for i in range(1,len(exon),2):
                rect = Rectangle((exon[i],trpNum-0.1),exon[i+1]-exon[i],0.2, color=col, fill=True)
                ax.add_patch(rect)
                if i < len(exon)-2:
                        ax.plot([exon[i+1],exon[i+2]],[trpNum,trpNum],color='black',linewidth=1)

trpNum = 0
for trpt,utr in GRCh38UTR['DDX11L1'].items():
        trpNum += 1
        if utr == []:
                continue
        for j in range(0,len(utr),2):
                rect = Rectangle((utr[j],trpNum-0.1),utr[j+1]-utr[j],0.2, color='black', fill=True)
                ax.add_patch(rect)
ax.set_xlabel('DDX11L1',color='blue',fontsize=14,fontweight='bold')
ax.yaxis.set_major_locator(ticker.FixedLocator(range(1,trpNum+1)))
ax.set_yticklabels(GRCh38Exon['DDX11L1'].keys())

plt.xlim(geneCord['DDX11L1'])
plt.ylim(0.5,trpNum+0.5)

plt.show()



figure_2.png
回复 支持 2 反对 0

使用道具 举报

0

主题

6

帖子

291

积分

中级会员

Rank: 3Rank: 3

积分
291
发表于 2017-2-18 20:04:24 | 显示全部楼层
本帖最后由 雨林木风 于 2017-2-25 20:13 编辑

思路:
1 构建有序字典geneCord记录gene的起至坐标,构建有序字典GRCh38Exon记录gene每个转录本的外显子的链信息和起止坐标,构建有序字典GRCh38UTR记录每个转录本的UTR起止坐标;
2 利用python的matplotlib包绘图

[Python] 纯文本查看 复制代码
# -*-coding:utf-8-*-

import re
import sys
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib.patches as patches
from matplotlib.patches import Rectangle
from collections import OrderedDict
matplotlib.style.use("ggplot")

args = sys.argv

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

def main():
    with open(args[1],'r') as fp_gtf:
        for line in fp_gtf:
            if line.startswith("#"):
                continue
            line_list = line.strip("\n").split("\t")
            chr = line_list[0]
            type = line_list[2]
            start = int(line_list[3])
            end = int(line_list[4])
            orientation = line_list[6]
            attr = line_list[8]

            if type == "gene":
                genename = re.search(r'gene_name "([^;]+)";?', attr).group(1)
                geneCord[genename] = [start,end]
                GRCh38Exon[genename] = OrderedDict()
                GRCh38UTR[genename] = OrderedDict()

            elif type == "transcript":
                trptname = re.search(r'transcript_name "([^;]+)";?', attr).group(1)
                genename = re.search(r'gene_name "([^;]+)";?', attr).group(1)
                if genename not in geneCord:
                    continue
                GRCh38Exon[genename][trptname] = []
                GRCh38Exon[genename][trptname].append(orientation)
                GRCh38UTR[genename][trptname] = []

            elif type == "exon":
                trptname = re.search(r'transcript_name "([^;]+)";?', attr).group(1)
                genename = re.search(r'gene_name "([^;]+)";?', attr).group(1)
                if genename not in geneCord:
                    continue
                if trptname not in GRCh38Exon[genename]:
                    continue
                GRCh38Exon[genename][trptname].extend([start,end])

            elif type in ["five_prime_utr","three_prime_utr"]:
                trptname = re.search(r'transcript_name "([^;]+)";?', attr).group(1)
                genename = re.search(r'gene_name "([^;]+)";?', attr).group(1)
                if genename not in geneCord:
                    continue
                if trptname not in GRCh38UTR[genename]:
                    continue
                GRCh38UTR[genename][trptname].extend([start,end])

    gene_model("SAMD11")

#绘制基因模型
def gene_model(genename):
    fig = plt.figure()
    ax = fig.add_subplot(111) #将画布分割成1行1列,图像画在从左到右从上到下的第1块

    trptNum = 0
    for trpt, exon in GRCh38Exon[genename].items():
        trptNum += 1
        if exon[0] == '+':
            col = 'limegreen'
        elif exon[0] == '-':
            col = 'magenta'
        for i in range(1,len(exon),2):
            #画外显子 [Rectangle((x,y), width, height)]
            rect = Rectangle((exon[i], trptNum-0.1), exon[i+1]-exon[i], 0.2, color=col, fill=True)
            ax.add_patch(rect)
            #添加内含子线
            if i < len(exon)-2:
                ax.plot([exon[i+1],exon[i+2]],[trptNum,trptNum],color='black',linewidth=1)

    trptNum = 0
    for trpt, utr in GRCh38UTR[genename].items():
        trptNum += 1
        if utr == []:
            continue
        for j in range(0,len(utr),2):
        #画UTR区 [Rectangle((x,y),width,height)]
            rect = Rectangle((utr[j], trptNum-0.1), utr[j+1]-utr[j], 0.2, color='black', fill=True)
            ax.add_patch(rect)

    ax.set_xlabel(genename, color='blue', fontsize=14, fontweight='bold')
    ax.yaxis.set_major_locator(ticker.FixedLocator(range(1,trptNum+1)))
    ax.set_yticklabels(GRCh38Exon[genename].keys(),fontweight='bold')

    plt.xlim(geneCord[genename])
    plt.ylim([0.5,trptNum+0.5])
    plt.show()

if __name__ == "__main__":
    main()


结果: figure_1.png
小结:python绘图还不是很熟悉,是看了老师的视频后仿写的,今后对于自己不熟悉的内容后续要加强提高google搜索学习的能力。
回复 支持 1 反对 0

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-2-18 22:28:17 | 显示全部楼层
034-perl+shell

数据为gencode.v7.annotation_goodContig.gtf

先用perl处理下数据,提取出关键数据即可,如gene_name为ANXA1的gene长度,以及其exon信息,perl代码如下:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my $in = "gencode.v7.annotation_goodContig.gtf";
my $out = "transcript.txt";
my $num = 0;

open(my $fh, $in) or die "Cannot open $in";
open(my $fh_out, ">$out") or die "Cannot write to $out";

while(<$fh>){
	chomp;
	next if ($_ =~ /^\#\#/);
	next if ($_ =~ /^\D/);
	my @array = split /\t/, $_;
	if($array[0] == 9 && $array[1] eq "HAVANA" && $array[8]=~/ANXA1/){
		if ($array[2] eq "gene"){
			print $fh_out "chr9\t";
		}elsif($array[2] eq "transcript"){
			$num++;
			print $fh_out "$num\t";
		}elsif($array[2] eq "exon"){
			print $fh_out "$num\t";
		}else{
			next;
		}
		print $fh_out "$array[2]\t$array[3]\t$array[4]\n";
	}
}
close $fh;
close $fh_out;


然后用R进行数据的最终处理,如去掉一些画图用不到的信息以及添加上颜色信息,最后用ggplot2进行画图,R代码如下:
[AppleScript] 纯文本查看 复制代码
library(ggplot2)

data <- read.table(file = "transcript.txt", sep = "\t", stringsAsFactors = F)

xlab_name <- paste(data[1,1], paste(data[1,3],data[1,4],sep = "-"), sep = ":")
title <- "ANXA1"
data <- data[-1,]
data <- data[data[,2]!="transcript",]

num <- length(unique(as.factor(data[,1])))
color <- rainbow(num)

data_min <- min(c(data[,3], data[,4]))
data$new_start <- data[,3] - data_min
data$new_end <- data[,4] - data_min
data$new_color <- color[as.numeric(data[,1])]

p <- ggplot()+geom_line()+
  annotate("segment", x = data$new_start, xend = data$new_end, y = data[,1], yend = data[,1], color = data$new_color,size = 2)+
  labs(x = xlab_name, y = "number", title = title)
p


图片如下:


Rplot.png
回复 支持 2 反对 0

使用道具 举报

0

主题

4

帖子

54

积分

注册会员

Rank: 2

积分
54
发表于 2017-2-19 17:09:52 | 显示全部楼层
x2yline 发表于 2017-2-15 19:37
77-R-python-shell x2yline

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

大神你的python代码没看太懂,就是为何要引入buffer?为何要写成
[AppleScript] 纯文本查看 复制代码
buffer = buffer[-1] + f.read(chunk)
print((buffer[buffer.find('\n')+1:buffer.find('\n')+3])

这样写的好处是什么?我感觉有点绕。
我只会最简单粗暴的for 循环迭代文件内容,然后跳过不需要的行
回复 支持 反对

使用道具 举报

2

主题

34

帖子

777

积分

高级会员

Rank: 4

积分
777
发表于 2017-2-19 20:56:16 | 显示全部楼层
zeyuchen2016 发表于 2017-2-19 17:09
大神你的python代码没看太懂,就是为何要引入buffer?为何要写成
[mw_shl_code=applescript,true]buffer  ...

看第二题的嘉宾视频就懂了,buffer读取比逐行快的不只是一点点。
回复 支持 反对

使用道具 举报

4

主题

24

帖子

184

积分

注册会员

Rank: 2

积分
184
发表于 2017-2-22 01:36:18 | 显示全部楼层
本帖最后由 tsznxx 于 2017-2-22 01:52 编辑

不到20行的Python 代码实现:
重要的事情说三遍:首先把GTF转成genepred格式(或bed12)!首先把GTF转成genepred格式(或bed12)!首先把GTF转成genepred格式(或bed12)!
[Python] 纯文本查看 复制代码
import matplotlib.pyplot as pltcolors = ['red','blue','purple']  # color loops
gids = []
for cnt,line in enumerate(open("/scratch/bcb/ywang52/TData/test/test.genepred")):
    gene = line.rstrip().split("\t")
    gids.append(gene[0])
    plt.plot([int(gene[3]),int(gene[4])],[cnt,cnt],linewidth=1,color='black') # draw intron
    txstart, txstop = int(gene[5]),int(gene[6])
    for start,stop in zip([int(s) for s in gene[8].rstrip(',').split(',')],[int(s) for s in gene[9].rstrip(',').split(',')]):
        if start > txstop or stop < txstart: # UTR: no overlap with ORF
            plt.plot([start,stop],[cnt,cnt],linewidth=3,color='grey')
        else: # draw ORF
            plt.plot([start,stop],[cnt,cnt],linewidth=3,color='grey')
            plt.plot([max(txstart,start),min(txstop,stop)],[cnt,cnt],linewidth=5,color=colors[cnt%len(colors)])                
plt.yticks(range(cnt+1),gids)
plt.ylim(-0.5,cnt+0.5)
plt.tight_layout()
plt.savefig("test.pdf")


结果如下图:
test.png



回复 支持 1 反对 0

使用道具 举报

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-2-23 15:19:14 | 显示全部楼层
tsznxx 发表于 2017-2-22 01:36
不到20行的Python 代码实现:
重要的事情说三遍:首先把GTF转成genepred格式(或bed12)!首先把GTF转成gene ...

厉害!
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 1 反对 0

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-20 21:07 , Processed in 0.034628 second(s), 22 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.