搜索
查看: 3682|回复: 39

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

[复制链接]

631

主题

1158

帖子

3885

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3885
发表于 2017-1-4 23:43:32 | 显示全部楼层 |阅读模式
转载自我自己的博客:画基因结构图!
可以下载各种gtf,从NCBI,ENSEMBL,UCSC,GENCODE都可以!(记住,你下载什么样的gtf就需要修改成什么样的代码!!!)
重点就是得到所有基因的转录本个数,以及每个转录本的外显子的坐标。
如图:


比如对这个ANXA1基因来说,非常多的转录本,但是基因的起始终止坐标,是所有转录本起始终止坐标的极大值和极小值!同时,它是一个闭合基因,因为它存在一个转录本的起始终止坐标等于该基因的起始终止坐标。可以看到它的外显子并不是非常整齐的,虽然多个转录本会共享某些外显子,但是也存在某些外显子比同区域其它外显子长的现象!

再比如下面这个例子:对DDX11L11这个基因来说,前两个外显子都不会翻译,直到第三个外显子才开始翻译,构成CDS。
有些转录本是没有utr的,所以该转录本的起始坐标,就是CDS的起始坐标
这个非常有用,可以更新自己的一些概念:
1,如果基因有多个转录本,基因的起始坐标,就是该基因所有转录本的第一个外显子的起始坐标的最小值,同理基因的终止坐标,就是该基因的所有转录本的最后一个外显子的终止坐标的最大值。
2,通过这个概念,可以把基因分成闭合基因和非闭合基因。 闭合基因:有一个最长转录本使得基因起始终止坐标等于该最长转录本的起始终止坐标。(这个是我乱说的,并没有这个定义)
3,如果基因只有一个转录本,那么基因的起始终止坐标,就是转录本的起始终止坐标!
4,一个基因的一个转录本的5’utr区域可以包括多个外显子区域,前者是翻译行为,后者是转录行为
5,起始密码子和终止密码子是CDS的起止处,是基于翻译的概念
6,一个基因的多个转录本的外显子坐标不一定会排列整齐,每个转录本的剪切位点并不一定要比其它转录本一致!



R实现的代码如下:
[AppleScript] 纯文本查看 复制代码
rm(list=ls())
## [url=http://www.broadinstitute.org/cancer/cga/sites/default/files/data/tools/rnaseqc/gencode.v7.annotation_goodContig.gtf.gz]http://www.broadinstitute.org/ca ... n_goodContig.gtf.gz[/url]
setwd('tmp')
gtf <- read.table('gencode.v7.annotation_goodContig.gtf.gz',stringsAsFactors = F,
                  header = F,comment.char = "#",sep = '\t'
                  )
table(gtf[,2])
gtf <- gtf[gtf[,2] =='HAVANA',]
gtf <- gtf[grepl('protein_coding',gtf[,9]),]

lapply(gtf[1:10,9], function(x){
  y=strsplit(x,';')  
})

gtf$gene <- lapply(gtf[,9], function(x){
  y <- strsplit(x,';')[[1]][5]
  strsplit(y,'\\s')[[1]][3]
  }
)

draw_gene = 'TP53'
structure = gtf[gtf$gene==draw_gene,]
 
colnames(structure) =c(
  'chr','db','record','start','end','tmp1','tmp2','tmp3','tmp4','gene'
)


gene_start <- min(c(structure$start,structure$end))
gene_end <- max(c(structure$start,structure$end))
tmp_min=min(c(structure$start,structure$end))
structure$new_start=structure$start-tmp_min
structure$new_end=structure$end-tmp_min
tmp_max=max(c(structure$new_start,structure$new_end))
num_transcripts=nrow(structure[structure$record=='transcript',])
tmp_color=rainbow(num_transcripts)

x=1:tmp_max;y=rep(num_transcripts,length(x))
#x=10000:17000;y=rep(num_transcripts,length(x))
plot(x,y,type = 'n',xlab='',ylab = '',ylim = c(0,num_transcripts+1))
title(main = draw_gene,sub = paste("chr",structure$chr,":",gene_start,"-",gene_end,sep=""))
j=0;
tmp_legend=c()
for (i in 1:nrow(structure)){
  tmp=structure[i,]
  if(tmp$record == 'transcript'){
    j=j+1
    tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
  }
  if(tmp$record == 'exon') lines(c(tmp$new_start,tmp$new_end),c(j,j),col=tmp_color[j],lwd=4)
}



我看了python讲师录制的视频,实在是太赞了,我都想赶快去学一下python了,恭喜我们的学员,有福了!
他的代码量比较多,暂时我们不建议把源码放在这里,你们根据视频慢慢敲吧,记住一步步的实现,自己动手,学的更扎实!







本帖子中包含更多资源

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

x



上一篇:外显子重复的基因,是如何定量的呢??
下一篇:生信编程直播第二题-hg19基因组序列的一些探究
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

1

主题

31

帖子

414

积分

中级会员

Rank: 3Rank: 3

积分
414
发表于 2017-2-15 19:37:03 | 显示全部楼层
本帖最后由 x2yline 于 2017-2-26 12:12 编辑

77-R-python-shell x2yline

数据来源:(第三题的文件)ftp://ftp.ensembl.org/pub/release-87/gtf/homo_sapiens/Homo_sapiens.GRCh38.87.chr.gtf.gz

1. R语言实现
[AppleScript] 纯文本查看 复制代码
get_gene_name_vect<- function(gtf_data){
  gene_info_vect <- gtf_data[,9]
  gene_info_vect[gtf_data[,3]!='gene'] = ''
  gene_name_vect <- sapply(gene_info_vect,function(x){
    strsplit(strsplit(x,';')[[1]][3], '\"')[[1]][2]
  })
  names(gene_name_vect) = NULL
  return(gene_name_vect)
}

find_target_data <- function(gene, all_data, gene_name_vect){
  target_row <- which(gene_name_vect == toupper(gene))
  i <- 1
  while( is.na(gene_name_vect[target_row +i])){
    i <- i+ 1
  }
  target_data<-all_data[target_row:(target_row+i-1),]
  colnames(target_data) =c('chr','db','record','start','end','tmp1','strand','tmp3','tmp4')
  return(target_data)
}

plot_gene_structure <- function (target_data, gene, rect_width=0.8, intron_width=0.2) {
  gene_start <- target_data$start[1]
  gene_end <- target_data$end[1]
  transcript_num <- length(which(target_data$record == 'transcript'))
  tmp_colors <- c('green', 'red', 'black', 'blue', 'yellow', 'blue','grey','grey')
  names(tmp_colors) <- c('gene', 'CDS', 'transcript', 'exon','start_codon','stop_codon','three_prime_utr','five_prime_utr')
  rect_sep <- c(intron_width,rep(rect_width,7))
  names(rect_sep) <-c('transcript','gene', 'exon', 'CDS','start_codon','stop_codon','three_prime_utr','five_prime_utr')
  par(mar=c(6,6,2,6), bty='n', new=F)
  plot(gene_start:gene_start, 1, type = 'n', xlab='', ylab ='',ylim = c(0,transcript_num*2+1), xlim = c(gene_start,gene_end), yaxt='n',xaxt='n')
  title(main = gene,sub = paste("chr",target_data$chr,": ",gene_start,"-",gene_end,sep=""))
  if (target_data$strand[1]=='-'){
    strand <- 1
  }else { strand <-2 }
  j <- 0
  y_lab <- c(j)
  rect_ybottom <- j-rect_sep[1]
  rect_ytop <- j+rect_sep[1]
  #rect(gene_start, rect_ybottom, gene_end, rect_ytop, col = tmp_colors['gene'], border = F)
  arrows(gene_start, j+rect_sep[1], gene_end, j+rect_sep[1], code= strand, col = tmp_colors['gene'], lwd = 5, angle=30)
  for (rows in (2:dim(target_data)[1])){
    if (target_data$record[rows] == 'transcript') {
      j <- j + rect_sep[2] * 2.5
      y_lab <- c(y_lab,j)
    }
    rect_ybottom <- j-rect_sep[target_data$record[rows]]
    rect_ytop <- j+rect_sep[target_data$record[rows]]
    #print(tmp_colors[target_data$record[rows]])
    rect(target_data$start[rows], rect_ybottom, target_data$end[rows], rect_ytop, col=tmp_colors[target_data$record[rows]], border = NA)
  }

  axis(1,c(gene_start,gene_start+(gene_end-gene_start)/5,gene_start+(gene_end-gene_start)*2/5,gene_start+(gene_end-gene_start)*3/5,gene_start+(gene_end-gene_start)*4/5,gene_end))

  axis(2,y_lab, labels=c('gene',paste('transcipt', seq(1:transcript_num), sep=' ')),las=1)
  legend(x=gene_end, transcript_num*(rect_sep[2] * 2.5)+1,horiz=F,
         box.lty=0,xpd=T,
         c('gene','non_CDS_exon','CDS_exon','UTR_exon','intron'),
         fill=c('green','blue','red','grey','black'),
         border = F,bty='n',
         cex=0.8)
}


setwd('E:\\r\\biotrainee_demo\\class5')
file <- 'Homo_sapiens.GRCh38.87.chr.gtf'

library(data.table)
#比较耗时
all_data<-as.data.frame(fread(file, header=F, sep='\t', stringsAsFactors=F, skip =5))
#load('all.Rdata')

gene='tp53'
gene='ANxa1'

#比较耗时
gene_name_vect <- get_gene_name_vect(all_data)

target_data <- find_target_data (gene, all_data, gene_name_vect)

png(paste(gene,'.png'),width = 1200, height = 680)
plot_gene_structure(target_data, gene)
dev.off()






(2)python3 实现:

更新后的代码如下:1. 加上了solid_capstyle = 'butt'这个参数,去掉了直线的‘帽子’

2.设置antialiased=False, 取消直线边缘模糊化处理

3.每段序列左右两端加了0.5
[Python] 纯文本查看 复制代码
def find_target_data(gene_name, gtf_file_path, chunk=3072*2048):
    with open(gtf_file_path) as f:
        header = ['chr','db','record','start','end','tmp','strand','tmp','info']
        target_data = {} 
        get = 0
        buffer_list =['']
        print('Please wait for 10 seconds: ')
        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')
                        for i in range(len(line_list)):
                            try:
                                target_data[header[i]].append(line_list[i])
                            except:
                                target_data[header[i]] = [line_list[i]]
            else:
                if get == 1:
                    break
            if len(buffer) < chunk:
                break
        if not target_data:
            print('\n\n There is some wrong with your gene name!\n')
            raise NameError('your gene_name is not exit')
        print ("\nHave got the gene information!")
    return(target_data)

def draw_gene_structure(gene_name,target_data, png_path='',line_width=5):
    gene_symbol = gene_name.upper()
    if not png_path:
        png_path = gene_symbol+'.png'
    #定义颜色的字典
    tmp_colors = ['lime','red', 'blue', 'yellow','yellow','w']
    names_tmp_colors = [ 'gene','CDS','exon','three_prime_utr','five_prime_utr','stop_codon']
    colors_legend_name = ['gene','CDS_exon', 'non_CDS_exon', 'UTR_exon']
    color_dict = dict(zip(names_tmp_colors,tmp_colors))

    # 提取转录本名称
    import re
    transcript_list = []
    for i in target_data['info']:
        try:
            transcript_name = re.findall('transcript_name "(.*?)"',i)[0]
            if transcript_name not in transcript_list:
                transcript_list.append(transcript_name)
        except:
            pass
    # 计算转录本数目
    transcript_num = 0
    for i in target_data['record']:
        if i =='transcript':
           transcript_num += 1

    import numpy as np
    import matplotlib.patches as mpatches
    import matplotlib.pyplot as plt
    import matplotlib.lines as lines
    fig = plt.figure(1)
    # 设置不透明度,默认为1
    fig.patch.set_alpha(1)
    fig.patch.set_facecolor('w')
    num = 0#当前转录本数目标志
    warnings = []
    for i in range(len(target_data['record'])):
        if target_data['record'][i] == 'gene':
            #判断正反链
            if target_data['strand'][i] == '+':
                arr = '->'
            else:
                arr = '<-'
            #图的第一个区域
            # add_axes 是在一张图上指定特定区域作图,第一个数字为从左边%20处,下面20%处开始,宽50%,高60%区域作图
            ax = fig.add_axes([0.2,0.2,0.5,0.6])
            # 定义基因方向箭头
            arrow= mpatches.FancyArrowPatch(
            (int(target_data['start'][i]), 0.1),
            (int(target_data['end'][i]), 0.1),
            arrowstyle= arr,
            mutation_scale=25, lw=1, color='lime',antialiased=True)#antialiased默认为True,边缘平滑处理
            # 画箭头
            ax.add_patch(arrow)
            # 坐标轴标签
            ax.set_xlim([int(target_data['start'][i]), int(target_data['end'][i])])
            ax.set_ylim([-0.5, transcript_num+1])
            ax.set_xticks(np.linspace(int(target_data['start'][i]),int(target_data['end'][i]),5))
            ax.set_yticks([0.1]+list(range(1,transcript_num +1)))
            ax.set_yticklabels(['gene']+transcript_list)
            ax.set_xticklabels([int(i) for i in np.linspace(int(target_data['start'][i]),int(target_data['end'][i]),5)])
            # 坐标轴显示
            ax.spines['top'].set_visible(False)
            ax.spines['left'].set_visible(False)
            ax.spines['right'].set_visible(False)
            ax.get_xaxis().tick_bottom()
            ax.get_yaxis().tick_left()
            ax.get_xaxis().set_tick_params(direction='out')
            ax.tick_params(axis=u'y', which=u'both',length=0)
            # 坐标轴字体大小
            for tick in ax.xaxis.get_major_ticks():
                tick.label.set_fontsize(6) 
            for tick in ax.yaxis.get_major_ticks():
                tick.label.set_fontsize(6) 
        elif target_data['record'][i] == 'transcript':
            num += 1 # 转录本所有区域计数作图
            line1 = [(int(target_data['start'][i]),num), (int(target_data['end'][i]),num)]
            (line1_xs, line1_ys) = zip(*line1)
            ax.add_line(lines.Line2D(line1_xs, line1_ys, linewidth=0.2,
                                     solid_capstyle = 'butt',solid_joinstyle='miter',
                                     antialiased=False,color='black'))
        elif target_data['record'][i] in color_dict.keys():
            # 添加结构图
            line2 = [(int(target_data['start'][i])-0.5,num), (int(target_data['end'][i])+0.5,num)]
            (line2_xs, line2_ys) = zip(*line2)
            ax.add_line(lines.Line2D(line2_xs, line2_ys, 
                                solid_capstyle = 'butt',solid_joinstyle='miter',
                                linewidth=int(line_width), alpha = 1,
                                color=color_dict[target_data['record'][i]],
                                antialiased=False))
        else:
            warnings.append(target_data['record'][i])
    if warnings:
        print('\nTips: ')
        print(' and '.join([i for i in set(warnings)]) + ' is not in our consideration!!!!!!')


    # 做图例
    # add_axes 是在一张图上指定特定区域作图,第一个数字为从左边%74处,下面20%处开始,宽20%,高60%区域作图
    ax_legend = fig.add_axes([0.76,0.2,0.2,0.6])
    #ax_legend.set_xticks([])
    #ax_legend.set_yticks([])
    for i in range(len(colors_legend_name)):
        line3 = [(0, (9-i)*0.1),(0.1, (9-i)*0.1)]
        (line3_xs, line3_ys) = zip(*line3)
        ax_legend.add_line(lines.Line2D(line3_xs, line3_ys, linewidth=5,
                                color=color_dict[names_tmp_colors[i]],
                                solid_capstyle = 'butt',solid_joinstyle='miter',
                                antialiased=False))
        ax_legend.text(0.2, (8.9-i)*0.1,colors_legend_name[i] , fontsize=6)
    ax_legend.set_axis_off()
    # 加标题
    fig.suptitle('\n\n\nchr' + str(target_data['chr'][0])+ ': ' + gene_symbol, fontsize=10)
    # 保存图片
    fig.savefig(png_path, dpi=150)
    plt.show()
    print('\nThe picture file is completed: ' + png_path)
    print("All transcripts of " + gene_name + ':\n' + " ".join(sorted(transcript_list)))
    return(png_path)



# gtf文件下载地址:
#[url]ftp://ftp.ensembl.org/pub/release-87/gtf/homo_sapiens/Homo_sapiens.GRCh38.87.chr.gtf.gz[/url]
gtf_file_path = 'Homo_sapiens.GRCh38.87.chr.gtf'

gene_name = input('Please enter the gene_name: ').strip()
line_width = input('Please enter the line_width: ').strip()

if not line_width:
    line_width = 5
if not gene_name:
	gene_name = input('Please enter the gene_name: ').strip()
else:
    gene_name = 'TP53'
    #gene_name = 'hoxc8'
    #gene_name = 'AnXA1'


#计算消耗时间
import time
start = time.time()

# 提取目标基因的数据
target_data = find_target_data(gene_name=gene_name, gtf_file_path=gtf_file_path)
# 作图
draw_gene_structure(gene_name,target_data, line_width=line_width)

print('\n used %.2f s to get the target information and get its structure!'%(time.time()-start))


新代码产生的图如下







总结:R语言读取大文件的速度太慢,不知道还有没有什么优化的方法;

python的read(buffer)速度非常赞,作图功能也是十分强大。




问题:
All transcripts of TP53:
TP53-001 TP53-002 TP53-003 TP53-004 TP53-005 TP53-006 TP53-007 TP53-008 TP53-009 TP53-013 TP53-014 TP53-015 TP53-016 TP53-017 TP53-018 TP53-019 TP53-020 TP53-021 TP53-022 TP53-023 TP53-024 TP53-025 TP53-026 TP53-027 TP53-028 TP53-029 TP53-201 TP53-202

发现没有TP53-010、TP53-011和TP53-012
查找验证发现:
http://asia.ensembl.org/Homo_sapiens/Gene/Summary?g=ENSG00000141510;r=17:7661779-7687550
中也没有名称为这三个的转录本









本帖子中包含更多资源

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

x
回复 支持 5 反对 0

使用道具 举报

95

主题

119

帖子

618

积分

版主

Rank: 7Rank: 7Rank: 7

积分
618
发表于 2017-2-26 13:43:29 | 显示全部楼层
002 +python +shell

这题对我来说太难,一步一步跟着   “77-R-python-shell x2yline” 学的
先把过程实现,然后慢慢消化





本帖子中包含更多资源

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

x
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 2 反对 0

使用道具 举报

25

主题

100

帖子

758

积分

高级会员

Rank: 4

积分
758
发表于 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


图片如下:


本帖子中包含更多资源

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

x
回复 支持 2 反对 0

使用道具 举报

6

主题

23

帖子

177

积分

注册会员

Rank: 2

积分
177
发表于 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()



本帖子中包含更多资源

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

x
回复 支持 2 反对 0

使用道具 举报

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

使用道具 举报

3

主题

21

帖子

128

积分

注册会员

Rank: 2

积分
128
发表于 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")


结果如下图:




本帖子中包含更多资源

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

x
回复 支持 1 反对 0

使用道具 举报

631

主题

1158

帖子

3885

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3885
 楼主| 发表于 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

使用道具 举报

0

主题

11

帖子

71

积分

注册会员

Rank: 2

积分
71
发表于 2017-1-13 21:25:04 | 显示全部楼层
感谢群主的代码,我把这一题重复了一遍,代码太经典,我没法改,所以就不上传了。
有个问题不明白
tmp_legend表示的每条转录本的信息
if(tmp$record == 'transcript'){
    j=j+1
    tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
  }
问题是我在画图中没有发现tmp_legend在图中应该位于什么地方。标题,副标题都是明确的。所以,不知道为啥要写这个代码?
回复 支持 1 反对 0

使用道具 举报

631

主题

1158

帖子

3885

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3885
 楼主| 发表于 2017-2-1 09:10:54 | 显示全部楼层
13235695523 发表于 2017-1-13 21:25
感谢群主的代码,我把这一题重复了一遍,代码太经典,我没法改,所以就不上传了。
有个问题不明白
tmp_lege ...

聪明,这个代码是我博客历史遗留问题,你去看我博客就明白了!
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

13

帖子

87

积分

注册会员

Rank: 2

积分
87
发表于 2017-2-7 22:19:21 | 显示全部楼层
本帖最后由 tingting 于 2017-2-7 22:22 编辑

121-R/python-婷

本来想在视频出来前先自己思考一下该如何做,结果发现是画图题,完全不会,还是坐等视频吧。

看到python都能画出这么漂亮的图,真想赶紧看视频啊!

不过趁着视频还没出来,赶紧把前四题视频上面的代码再默写一遍巩固好,这样后面的题目才会越做越顺。

回复 支持 反对

使用道具 举报

631

主题

1158

帖子

3885

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
3885
 楼主| 发表于 2017-2-9 23:20:46 | 显示全部楼层
tingting 发表于 2017-2-7 22:19
121-R/python-婷

本来想在视频出来前先自己思考一下该如何做,结果发现是画图题,完全不会,还是坐等视频 ...

你的想法很赞,这是一个系列,前面的知识的顺利掌握是很重要的!加油
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

4

帖子

54

积分

注册会员

Rank: 2

积分
54
发表于 2017-2-13 15:27:52 | 显示全部楼层
没有搞懂这句代码:
[AppleScript] 纯文本查看 复制代码
gtf <- gtf[gtf[,2] =='HAVANA',]

难道“ENSEMBL”中的不要了?
回复 支持 反对

使用道具 举报

0

主题

15

帖子

147

积分

注册会员

Rank: 2

积分
147
发表于 2017-2-14 14:30:19 | 显示全部楼层
弱弱的问一句,第五题的讲解是不是没有出完啊 ? 讲perl的那个视频只制作了数据提取部分的呢 ?
回复 支持 反对

使用道具 举报

0

主题

15

帖子

147

积分

注册会员

Rank: 2

积分
147
发表于 2017-2-14 16:27:52 | 显示全部楼层
zeyuchen2016 发表于 2017-2-13 15:27
没有搞懂这句代码:
[mw_shl_code=applescript,true]gtf

可能就是,我刚看了R的教学视频,里面提到过 说HAVANA比ENSEMBL要全一些,所以选择性的忽略掉了这个数据库的,就酱
回复 支持 反对

使用道具 举报

0

主题

4

帖子

54

积分

注册会员

Rank: 2

积分
54
发表于 2017-2-15 12:12:45 | 显示全部楼层
Aiyawq 发表于 2017-2-14 16:27
可能就是,我刚看了R的教学视频,里面提到过 说HAVANA比ENSEMBL要全一些,所以选择性的忽略掉了这个数据 ...

还没看视频呢!谢谢
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|关于我们|手机版|小黑屋|生信技能树    

GMT+8, 2017-10-18 11:39 , Processed in 0.046012 second(s), 39 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.