搜索
楼主: Jimmy

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

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

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

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





本帖子中包含更多资源

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

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

使用道具 举报

2

主题

21

帖子

241

积分

中级会员

Rank: 3Rank: 3

积分
241
发表于 2017-2-26 22:14:04 | 显示全部楼层
本帖最后由 AdaWon 于 2017-2-26 22:37 编辑

147- python
1. 理解了第一部分的函数,即厘清gene/transcript/exon/UTR 之间的层级关系,并利用有序字典,构建基因、对应的转录本、对应的外显子以及UTR的坐标和相应的联系;2. 学习matplotlib,但是运行不成功,总是报错,目前还未解决,暂搁置,MARK,以后解决!!!
I would appreciate it if someone could help me!!!
参考链接1  patches
参考链接2  Matplotlib_slides
参考链接3  matplotlib-user-guide

[Python] 纯文本查看 复制代码
####parse exon coordinates of each transcript
from collections import OrderedDict 
##导入特定的成员
#Python拥有一些内置的数据类型,比如str, int, list, tuple, dict等, collections 模块在这些内置数据类型的基础上,提供了几个额外的数据类型:
##namedtuple(): 生成可以使用名字来访问元素内容的tuple子类
##deque: 双端队列,可以快速的从另外一侧追加和推出对象
##Counter: 计数器,主要用来计数
##OrderedDict: 有序字典
##defaultdict: 带有默认值的字典

geneCord = OrderedDict()
GFF_Exon = OrderedDict()
GFF_UTR = OrderedDict()

with open('/data0/home/guixiang/sgn_genome/ITAG3.0_gene_models.gff','rt') as f:
        #r、w、a为打开文件的基本模式,对应着只读、只写、追加模式;
        #b、t、+、U这四个字符,与以上的文件打开模式组合使用,二进制模式,文本模式,读写模式、通用换行符
        for line in f:
                line = line.rstrip()  #rstrip() 删除 string 字符串末尾的指定字符(默认为空格)
                if line.startswith('#'):
                        continue

                lst = line.split('\t')

                if lst[-1].startswith('ID=gene'):
                        geneInfo = lst[-1].split(';')
                        geneName = geneInfo[0].split('=')[1]         #geneInfo[0] is ID=gene:Solyc00g005000.3;
                        geneCord[geneName] = [int(lst[3]),int(lst[4])] #start,end
                        GFF_Exon[geneName] = OrderedDict()
                        GFF_UTR[geneName] = OrderedDict()

                elif lst[-1].startswith('ID=mRNA'):
                        trptInfo = lst[-1].split(';')
                        trptName = trptInfo[0].split('=')[1]        #trptInfo[0] is ID=mRNA:Solyc00g005000.3.1;
                        GFF_Exon[geneName][trptName] = []        #字典嵌套
                        GFF_Exon[geneName][trptName].append(lst[6]) #strand
                        GFF_UTR[geneName][trptName] = []

                elif lst[2] == 'exon':
                        GFF_Exon[geneName][trptName].extend([int(lst[3]),int(lst[4])])

                elif lst[2]        in ['three_prime_UTR','five_prime_UTR']:
                        GFF_UTR[geneName][trptName].extend([int(lst[3]),int(lst[4])])


###Failed part!!!
[Python] 纯文本查看 复制代码
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')

def gene_model(geneName):
        #标准用法是创建一个Figure实例,使用Figure创建一个或多个Axes或Subplot实例,并使用Axes实例的辅助方法来创建基本类型。
        fig = plt.figure() #创建一个Figure实例
        ax = fig.add_subplot(111) #numrows,numcols,fignum;其中fignum的范围是从1到numrows * numcols


        trptNum = 0
        for trpt, exon in GFF_Exon[geneName].items():  
        #Python 字典(Dictionary) items() 函数以列表返回可遍历的(键, 值) 元组数组。
                trptNum += 1
                if exon[0] == '+':
                        col = 'limegreen'
                elif exon[0] == '-':
                        col = 'magenta'

                for i in range(1,len(exon),2):  
                        #len() 方法返回列表元素个数
                        #range(1,5,2) #代表从1到5,间隔2(不包含5)
                        #draw exons
                        rect = Rectangle((exon[i],trptNum-0.1),exon[i+1]-exon[i],0.2,color=col,fill=True)
                        #matplotlib.patches.Rectangle(xy, width, height, angle=0.0, **kwargs)
                        #Draw a rectangle with lower left at xy = (x, y) with specified width and height.
                        ax.add_patch(rect)

                        ##add intron lines
                        if i < len(exon)-2:
                                intronLength = exon[i+2]-exon[i+1]
                                #ax.plot([exon[i+1],exon[i+2]],[trptNum,trptNum],color = 'black',linewidth=1)


                                ### arrow line indicating the strand infomation
                                if exon[0] == '+':
                                        if intronLength < 500:
                                                # matplotlib.patches.Arrow(x, y, dx, dy, width=1.0, **kwargs)
                                                #Draws an arrow, starting at (x, y), direction and length given by (dx, dy) the width of the arrow is scaled by width;head_length:length of the arrow head
                                                ax.arrow(exon[i+1],trptNum,intronLength*0.9,0,head_width=0.03,head_length=intronLength*0.1,fc='gold',ec='gold')
                                        else:
                                                ax.arrow(exon[i+1],trptNum,intronLength-50,0,head_width=0.03,head_length=50,fc='gold',ec='gold')
                                elif exon[0] == '-':
                                        f intronLength < 500:
                                                ax.arrow(exon[i+2],trptNum, -intronLength*0.9,0,head_width=0.03,head_length=intronLength*0.1,fc='gold',ec='gold')
                                        else:
                                                ax.arrow(exon[i+2],trptNum, -intronLength+50,0,head_width=0.03,head_length=50,fc='gold',ec='gold')

        
        ### fill utr region with black
        trptNum = 0                                                                        
        for trpt,utr in GFF_UTR[geneName].items():
                trptNum += 1

                if utr == []:
                        continue

                for j in range(0,len(utr),2):
                        ## draw utrs [Rectangle((x,y),height,width)]
                        rect = Rectangle((utr[j],trptNum-0.1),utr[j+1]-utr[j],0.2, color='black',fill = Ture)
                        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(GFF_Exon[geneName].keys(),fontweight='bold')

        plt.xlim(geneCord[geneName])
        plt.ylim([1,trptNum+0.5])

        plt.show()





本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

0

主题

5

帖子

36

积分

新手上路

Rank: 1

积分
36
发表于 2017-2-27 09:42:25 | 显示全部楼层
这个,我能缓一缓吗 ,有点难的感觉
回复 支持 反对

使用道具 举报

0

主题

15

帖子

151

积分

注册会员

Rank: 2

积分
151
发表于 2017-2-27 11:08:32 | 显示全部楼层
本帖最后由 Aiyawq 于 2017-2-27 11:15 编辑
anlan 发表于 2017-2-18 22:28
034-perl+shell

数据为gencode.v7.annotation_goodContig.gtf

厉害了,我也是用perl先预处理的,但是一直弄不好后面R作图的那部分内容呢~
回复 支持 反对

使用道具 举报

0

主题

16

帖子

149

积分

注册会员

Rank: 2

积分
149
发表于 2017-2-28 12:28:11 | 显示全部楼层
本帖最后由 xxnku 于 2017-2-28 15:10 编辑

[Python] 纯文本查看 复制代码
#encoding=utf-8
from collections import OrderedDict
import sys
args = sys.argv

geneCord = OrderedDict()
GRch38Exon = OrderedDict()
GRch38UTR = OrderedDict()

with open(args[1]) as f:
    for line in f:
        line = line.rstrip()
        if line.startswith('#'):
            continue
        lst = line.split('\t')

        if lst[-1].startswith('ID=gene'):
            geneInfo = lst[-1].split(';')
            geneName = geneInfo[1].split('=')[1]
            geneCord[geneName] = [int(lst[3]),int(lst[4])]
            GRch38Exon[geneName] = OrderedDict()
            GRch38UTR[geneName] = OrderedDict()

        elif lst[-1].startswith('ID=transcript'):
            trptInfo = lst[-1].split(';')
            trptName = trptInfo[2].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])])
        
        elif lst[2] in ['three_prime_UTR','five_prime_UTR']:
            GRch38UTR[geneName][trptName].extend([int(lst[3]),int(lst[4])])

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()
#返回一个Figure对象,也就是一个图片

ax = fig.add_subplot(111)
#将画布分割成1行1列,图像画在从左到右从上到下的第一块(这里也就是画在中间)

trpNum = 0
for trpt,exon in GRch38Exon['SAMD11'].items():
    trpNum += 1
    if exon[0] == '+':
        col = 'blue'
    elif exon[0] == '-':
        col = 'red'
    for i in range(1,len(exon),2):
        #draw exons
        rect = Rectangle((exon[i],trpNum-0.1),exon[i+1]-exon[i],0.2,color=col,fill=True)
        ax.add_patch(rect)
        #Rectangle((x,y),width,height,fill=False/True:设置backgroud),  y设置成trpNum-0.1和height设置成0.2,
        #正好可以使exon和inton位于同一水平线上。

        #draw intron lines
        if i < len(exon)-2:
            x_values = [exon[i+1],exon[i+2]]
            y_values = [trpNum,trpNum]
            ax.plot(x_values, y_values, color = 'black',linewidth = 1)
            #ax.plot(x,y,color=,linewidth=,.....)

trpNum = 0
for trpt,utr in GRch38UTR['SAMD11'].items():
    trpNum += 1
    if not utr == []:
        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('The structures of SAMD11',color='blue',fontsize=14,fontweight='bold')
#设置x轴 label, fontweight文本粗细属性,normal:标准,bold:加粗,bolder:更粗,lighter:更细

ax.yaxis.set_major_locator(ticker.FixedLocator(range(1,trpNum+1)))
#设置y轴主刻度

ax.set_yticklabels(GRch38Exon['SAMD11'].keys())
#设置y轴 刻度label

plt.xlim(geneCord['SAMD11'])
#设置x轴刻度范围

plt.ylim(0.5,trpNum+0.5)
#设置y轴刻度范围

plt.show()
#打开matplotlib查看器,并显示绘制的图形

[i][i]

谢谢Python-东




本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

积分
559
QQ
发表于 2017-3-5 16:18:19 | 显示全部楼层
本帖最后由 旭日早升 于 2017-3-5 21:37 编辑

201-python-无名

这次作业拖太久了,实在是不应该,刚开始想做,却发现有点超出能力了,只能先试试,再看老师的视频学习了。这一讲以学习老师代码为主。
问题解析:问题说明:根据gtf文件注释信息画出自己感兴趣的基因的转录本结构。文件格式:gtf文件,详见 http://asia.ensembl.org/info/website/upload/gff.html,不同数据可能略有差异。问题解析:找到基因的所有转录本,以及每个转录本的外显子坐标信息。数据结构:GSTT1: {1: [0-1,1-2,2-3], 2: [0-1,2-3]}。提高效率:暂时未知。
这一讲由于不会画图,因而尝试目标定位在能把所有gene的转录本搞出来即可。由于gtf的结构一般现实gene,再跟着transcript,然后就是每个transcript的exon,所以我会每次判断一次type和name,再一次将他们存入字典。
我的尝试代码:
[Python] 纯文本查看 复制代码
import re, os
from collections import OrderedDict
from operator import itemgetter

os.chdir("./")

gene_transcript = OrderedDict()
with open("gencode.ucsc_v19.annotation.gtf",'rt') as f:
        i = 1
        for line in f:
                if i > 1000: #only test the first 1000 line
                        break
                i += 1 
                if line.startswith("#"):
                        continue
                line = line.rstrip()
                lst = line.split("\t")
                g_type = lst[2]
                g_start = lst[3]
                g_end = lst[4]
                g_name = lst[-1].split()[9]
                t_id = lst[-1].split()[3]
                if g_name not in gene_transcript and g_type == "gene":
                        gene_transcript[g_name] = {}
                        continue
                if g_name in gene_transcript and g_type == "transcript":
                        gene_transcript[g_name][t_id] = []
                        continue
                if g_name in gene_transcript and g_type == "exon" and t_id in gene_transcript[g_name]:
                        gene_transcript[g_name][t_id].append(str(g_start)+"-"+str(g_end))

print(gene_transcript)

学习了东老师的视频,东老师说的挺好的,他自己也不怎么用python画图,但是他可以利用已有知识和google的搜索功能,很快找到并学会了画图的工具包,这本身就是很好的能力,我结合老师代码对于我自己的数据和问题做了修改,一方面,我使用的是gencode的gtf文件,另一方面我是在服务器跑的数据,绘图是在我自己电脑上做的。
代码如下:
[Python] 纯文本查看 复制代码
#### Parse exon cordinates of each transcript
from collections import OrderedDict

geneCord = OrderedDict()
hg19Exon = OrderedDict()
hg19UTR = OrderedDict()

with open("gencode.ucsc_v19.annotation.gtf", 'rt') as f:
	for line in f:
		line = line.rstrip()
		if line.startswith("#"):
			continue

		lst = line.split("\t")
		if lst[2] == "gene":
			geneInfo = lst[-1].split(";")
			geneName = geneInfo[4].split("\"")[1]
			geneCord[geneName] = [int(lst[3]), int(lst[4])]
			hg19Exon[geneName] = OrderedDict()
			hg19UTR[geneName] = OrderedDict()
		elif lst[2] == "transcript":
			trptInfo = lst[-1].split(";")
			trptName = trptInfo[7].split("\"")[1]
			hg19Exon[geneName][trptName] = []
			hg19Exon[geneName][trptName].append(lst[6])
			hg19UTR[geneName][trptName] = []
		elif lst[2] == "exon":
			hg19Exon[geneName][trptName].extend([int(lst[3]), int(lst[4])])
		elif lst[2] == "UTR":
			hg19UTR[geneName][trptName].extend([int(lst[3]), int(lst[4])])

print(geneCord["GSTT1"])
print(hg19Exon["GSTT1"])
print(hg19UTR["GSTT1"])

#### output to draw a picture in my own computer
with open("geneCord.txt",'wt') as f:
	for gene, cord in geneCord.items():
		f.write(str(gene)+"\t"+"\t".join(map(str, cord))+"\n")
with open("hg19Exon.txt",'wt') as f:
	for gene, transcripts in hg19Exon.items():
		for transcript, cord in transcripts.items():
			f.write(str(gene)+"\t"+transcript+"\t"+"\t".join(map(str,cord))+"\n")
with open("hg19UTR.txt", 'wt') as f:
	for gene, UTRs in hg19UTR.items():
		for UTR, cord in UTRs.items():
			f.write(str(gene) + "\t"+UTR+"\t"+"\t".join(map(str,cord))+"\n")
#### read data in my own picture
geneCord = OrderedDict()
hg19Exon = OrderedDict()
hg19UTR = OrderedDict()

with open("geneCord.txt",'rt') as f:
	for line in f:
		line = line.rstrip()
		lst = line.split("\t")
		geneCord[lst[0]] = lst[1:]
with open("hg19Exon.txt",'rt') as f:
	for line in f:
		line = line.rstrip()
		lst = line.split("\t")
		if lst[0] not in hg19Exon:
			hg19Exon[lst[0]] = OrderedDict()
			hg19Exon[lst[0]][lst[1]] = lst[2:]
		else:
			hg19Exon[lst[0]][lst[1]] = lst[2:]

with open("hg19UTR.txt", 'rt') as f:
	for line in f:
		line = line.rstrip()
		lst = line.split("\t")
		if lst[0] not in hg19UTR:
			hg19UTR[lst[0]] = OrderedDict()
			hg19UTR[lst[0]][lst[1]] = lst[2:]
		else:
			hg19UTR[lst[0]][lst[1]] = lst[2:]
####picture (in python of my own computer)
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)

trptNum = 0
for trpt, exon in hg19Exon["GSTT1"].items():
	trptNum += 1 
	if exon[0] == "+":
		col = "limegreen"
	elif exon[0] == "-":
			col = "magenta"

	for i in range(1, len(exon), 2):
		## draw exons [Rectangle(x,y), height, width]
		rect = Rectangle((int(exon[i]), trptNum-0.1),int(exon[i+1])-int(exon[i]), 0.2, color=col, fill=True)
		ax.add_patch(rect)

		## add intro lines
		if i < len(exon)-2:
			ax.plot([int(exon[i+1]), int(exon[i+2])], [trptNum, trptNum], color="black", linewidth=1)

trptNum = 0 
for trpt, utr in hg19UTR["GSTT1"].items():
	trptNum += 1 
	if utr == []:
		continue
	for j in range(0, len(utr), 2):
		## draw UTRs [Rectangle((x,y), height, width)]
		rect = Rectangle((int(utr[j]), trptNum-0.1), int(utr[j+1])-int(utr[j]), 0.2, color="black", fill=True)
		ax.add_patch(rect)

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

plt.xlim(map(int,geneCord["GSTT1"]))
plt.ylim([0.5,trptNum+0.5])

plt.show()



感谢老师的视频录制,教会了我学习的态度和方法。




本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

10

主题

52

帖子

559

积分

版主

Rank: 7Rank: 7Rank: 7

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

厉害了我的哥,看起来这种文件格式需要处理的东西很少。
回复 支持 反对

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-3-6 16:35:50 | 显示全部楼层
R画图不是很懂,只能看老师的视频,然后学习下了!还有就是,图例太大,会覆盖掉一部分内容,暂时不知道怎么调整!
[AppleScript] 纯文本查看 复制代码
gtf <- read.table('gencode.v7.annotation_goodContig.gtf.gz',stringsAsFactors = F,
                  header = F,comment.char = "#",sep = '\t'
                  )##read.table读取压缩包文件,理解下各个参数
table(gtf[,2])
gtf <- gtf[gtf[,2] =='HAVANA',]##选择HAVANA数据库 当然也可以都保留(ensembl)
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]
  }
)##怎么把基因名字找出来 先用;分割,取第一个list第五个元素[[1]][5],在用空格分割,取第一个list第三个元素[[1]][3],处理完成,数据变成10列
 
draw_gene = 'TP53'
structure = gtf[gtf$gene==draw_gene,]##选择作图基因,这里是TP53
  
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',])##看有多少transcripts
tmp_color=rainbow(num_transcripts)##根据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=""))
	#legend('topleft',legend=tmp_legend,lty=1,lwd = 4,col = tmp_color)
  }
  if(tmp$record == 'exon') lines(c(tmp$new_start,tmp$new_end),c(j,j),col=tmp_color[j],lwd=4)
}
#legend('topleft',legend=tmp_legend,lty=1,lwd = 4,col = tmp_color)加了这行代码,画出的图例太大了,会覆盖掉一部分内容

看来jimmy的教程,又收获很多,受益匪浅,还需要继续努力!
027-R+Perl-游戏

本帖子中包含更多资源

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

x
人生若只如初见!
回复 支持 反对

使用道具 举报

5

主题

32

帖子

484

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
484
发表于 2017-3-8 09:45:05 | 显示全部楼层
先用perl处理文本
#! usr/bin/perl
use warnings;
use strict;
my ($gene_name,);
open INPUT,"gencode.v7.annotation_goodContig.gtf"; # 原始gtf文件;
# open TEST,">ANXA1_gene_structure.txt"; # 筛选数据;
open OUTPUT,">graph_data4.txt"; # 绘图数据;
$gene_name = 'ANXA1'; # 基因名;
my $num = 0;
while(<INPUT>){
    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 OUTPUT "chr9\t";
        }elsif($array[2] eq "transcript"){
            $num++;
            print OUTPUT "$num\t";
        }elsif($array[2] eq "exon"){
            print OUTPUT "$num\t";
        }else{
            next;
        }
        print OUTPUT "$array[2]\t$array[3]\t$array[4]\n";
    }
}
再用R绘图
data <- read.table("graph_data4.txt", sep = "\t", stringsAsFactors = F)
xlab_name <- paste(data[1,1], paste(data[1,3],data[1,4],sep = "-"), sep = ":")
data <- data[-1,]
num <- length(unique(as.factor(data[,1])))
color <- rainbow(num)
num <- length(unique(as.factor(data[,1])))
color <- rainbow(num)

data_min <- min(c(data[,3], data[,4]))
data_max<-max(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])]
x=seq(data_min,data_max,length=5)
y=rep(num,5)
plot(x,y,type="n",main = "ANXA1",xlab = xlab_name,ylim = c(0,num+1))
j=0;
for(i in 1:nrow(data)){
  tmp=data[i,]
  if(tmp$V2=="transcript")
  {
    j=j+1
    next
  }
  else
    {lines(c(tmp$V3,tmp$V4),c(j,j),lwd=4,col=color[j])}
  
}
回复 支持 反对

使用道具 举报

0

主题

6

帖子

89

积分

注册会员

Rank: 2

积分
89
发表于 2017-3-11 11:56:27 | 显示全部楼层
本帖最后由 banapi 于 2017-3-11 11:57 编辑

158-perl-shell-python
感叹一下perl处理文本好强大,我也更深刻的理解了gtf文件的格式,模仿代码:
[Perl] 纯文本查看 复制代码
#! urs/bin/perl
use warnings;
use strict;

my($row,$gene_name,$gene_start,$gene_end,$chr,@cells,$x,$count);

open INPUT,"gencode.v7.annotation_goodContig.gtf";
open OUTPUT,">shuchu.txt";

$gene_name = 'ANXA1';
while($row = <INPUT>){
        chomp($row);
        if(($row=~ m/HAVANA/)&&($row=~ m/gene_type "protein_coding";/)&&($row=~ m/gene_name "$gene_name";/)){
        @cells = split /\t/, $row;
        if($cells[2]=~ m/gene/){        
                $chr = $cells[0];
                $gene_start = $cells[3];
                $gene_end = $cells[4];
}
        elsif($cells[2]=~ m/transcript/){
                $count += 1;
}
        elsif($cells[2]=~ m/exon/){
                for($x=$cells[3]-$gene_start;$x<=$cells[4]-$gene_start;$x+=1){
                        print OUTPUT "$x\t$count\n";}
}

}
}
print "chr$chr\tgene $gene_name:\t$gene_start-$gene_end\n";

输出:
chr9        gene ANXA1:        75766673-75785309


本帖子中包含更多资源

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

x
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-16 15:44 , Processed in 0.049236 second(s), 22 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.