生信编程5. 根据GTF文件绘制多个转录本结构
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享：
生信菜鸟入门需要跨过的3个障碍：
了解数据背后的生物学意义，每个生物学数据都有背后的生物学意义 编写代码的能力 使用公共软件的能力
1. 使用R语言绘制
rm(list = ls())
gtf_1 <- read.table("/Users/duanshumeng/段舒檬/曾老师任务/生信编程/gencode.v32.annotation.gtf",sep = "\t",
stringsAsFactors = F,header = F)
gtf <- gtf_1[1:10000,]
#查看基因的来源，如果是多个数据库，取一个就够了
table(gtf[,2])
gtf <- gtf[gtf[,2]=="HAVANA",]
#只想要蛋白质编码的蛋白
gtf<- gtf[grepl("protein_coding",gtf[,9]),]
#取出基因的名称，关键步骤
library(stringr)
gtf$gene = lapply(gtf[,9], function(x){
if (length(unlist(strsplit(x,';')))==6){
str_split(str_split(x,'\\;',simplify = T)[,3],'\\ ',simplify = T)[,3]
} else{
str_split(str_split(x,'\\;',simplify = T)[,4],'\\ ',simplify = T)[,3]
}
})
#将绘图步骤包装成一个函数，可以多次使用
draw_gene_graph <- function(draw_gene,gtf){
structure = gtf[gtf$gene==draw_gene,]
colnames(structure) = c(
'chr','db','record','start','end','tmp1','tmp2','tmp3','tmp4','gene'
)
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))
plot(x,y,type = 'n',xlab = '',ylab = '',ylim = c(0,num_transcripts+1))
title(main = draw_gene,sub = paste(structure[structure$record=="gene",]$chr,":",
structure[structure$record=="gene",]$start,
"-",structure[structure$record=="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$end,sep = ""))
}
if(tmp$record == "exon") lines(c(tmp$new_start,tmp$new_end),c(j,j),
col=tmp_color[j],lwd=5)
}
}
#挑选一个基因绘图
draw_gene_graph("SAMD11",gtf)
结果图：
2. 使用python的gffutils绘制
gffutils介绍
gffutils是一个使用分层方式对GFF和GTF文件进行处理python模块。
pip3 install gffutils
pip3 install bcbio-gff
#1. 从gtf中提取目的基因的信息
grep -iw "DDX11L16" gencode.v32.annotation.gtf > DDX11L16.gtf
#2. 读取目标gtf文件
import gffutils
import os
import matplotlib.pyplot as plt
os.listdir()
fn = gffutils.example_filename(os.getcwd()+'/'+'DDX11L1.gtf')
db = gffutils.create_db(fn,dbfn='DDX11L1.db',force=True,keep_order=True,merge_strategy='merge',sort_attribute_values=True)
db = gffutils.FeatureDB('DDX11L1.db',keep_order=True)
gene = db['ENSG00000227159.8_PAR_Y'] #根据基因ID取出基因
gene.start
gene.end
transcripts = []
for i in db.children(gene,featuretype='transcript',order_by='start'):
transcripts.append(i)
exons = []
for i in db.children(gene,featuretype='exon',order_by='start'):
exons.append(i)
for i in exons:
plt.hlines(y=1,xmin=i.start-gene.start, i.end=gene.end)
plt.show()