生信编程3.hg38每条染色体的基因、转录本分布
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享，上一期是：生信编程系列（1-2）
首先我们需要的hg38的gtf文件
下载地址
wget ftp://ftp.ensembl.org/pub/release-87/gtf/homo_sapiens/Homo_sapiens.GRCh38.87.chr.gtf.gz
##gtf文件的简单介绍
第一列：表示该特征元素所在的染色体号，如果没有组装到染色体级别，那有可能是scaffold， 或者contig。这是基因组汇总的坐标系统，后续的一切注释信息都基于这一列。 第二列：表示这一行注释信息的来源 第三列：表示注释的类型，有基因，转录本，外显子，CDs等等 第四列：起始位置 第五列：终止位置 第六列：得分，改行特征的注视得分，如序列的相似性，如果没有对应的分数，一般使用.代替。 第七列：链，表示该行表示的特征所在的正负链，如果是.则表示不确定正负链，或者与链无关 第八列：相位，与蛋白质编码有关，一般用于CDS，值的范围从0-2，表示编码时阅读框的移动相位。 第九列：属性，对更多的注释信息进行描述
然后我们对gtf文件进行探索 ##首先是获得每条染色体含有的基因数量统计
第一种方法：使用perl进行统计
cat gencode.v30.annotation.gtf | perl -alne '{print if $F[2] eq "gene"}' | cut -f 1 | sort | uniq -c
==解释== $F[2]取文件的第三列，如果第三列的字符为“gene”的话，就输出该行，然后用cut以tab键作为分隔符进行分隔，取分隔后的第一个字符，也就是染色体号。最后对染色体出现的次数进行统计，就可以得到结果
第二种方法：使用awk进行统计
awk '{if($3=="gene") print $1}' gencode.v30.annotation.gtf | sort | uniq -c
得到的结果是一致的
查看基因的分类
cat gencode.v30.annotation.gtf | perl -alne '{next unless $F[2] eq "gene" ;/gene_type "(.*?)";/;print $1}' |sort |uniq -c
==解释== 这里比较重要的就是用到perl的正则表达/gene_type "(.*?)";/,括号里的是需要返回的内容
##使用python进行处理 首先对GTF的关键信息进行结构化
{"chr_id":{
"id":"",
"start":"value",
"end":"value",
"length":"value",
"gene":{
"geneid":"value",
"start":"value",
"end":"value",
"chr":"value",
"orientation":"value",
"length":"value",
"transcript":{
"transcript_id":"value",
"id":"value",
"start":"value",
"end":"value",
"chr":"value",
"parent":"value",
"length":"value",
"exon":[
{
"start":"value",
"end":"value",
"chr":"value",
"gene":"value",
"parent":"value",
"length":"value"
}],
"cds":[{
"start":"value",
"end":"value",
"chr":"value",
"gene':"value",
"parent":"value",
"length":"value"}]
}}}}
探究的内容：
染色体上基因数量的分布 基因长度的分布情况 基因的转录本数量分布 转录本的外显子数量统计 外显子的分布情况 探索的结果进行绘图
python代码实现
将整个过程分为4步
定义四个类，其中Genome_info为其他几个类的父类，其他几个类分别继承Genome_info的属性
#!/usr/bin/python
import sys
import re
args = sys.argv
class Genome_info:
def __init__(self):
self.chr = ""
self.start = 0
self.end = 0
class Gene(Genome_info):
def __init__(self):
Genome_info.__init__(self)
self.orientation = ""
self.id = ""
class Transcript(Genome_info):
def __init__(self):
Genome_info.__init__(self)
self.id = ""
self.parent = ""
class Exon(Genome_info):
def __init__(self):
Genome_info.__init__(self)
self.parent = ""
定义主函数
def main(args):
"""
input gtf file
"""
list_chr = []
list_gene = {}
list_transcript = {}
list_exon = []
with open(args[1]) as fp_gtf:
for line in fp_gtf:
if line.startswith("#"):
continue
lines = line.strip("\n").split("\t")
chr = lines[0]
type = lines[2]
start = int(lines[3])
end = int(lines[4])
orientation = lines[6]
attr = lines[8]
if not re.search(r'protein_coding',attr):
continue
if not chr in list_chr:
list_chr.append(chr)
if type == "gene":
gene = Gene()
id = re.search(r'gene_id "([^;]+)";?', attr).group(1)
gene.chr = chr
gene.start = start
gene.end = end
gene.id = id
gene.orientation = orientation
list_gene[id] = gene
print(id)
elif type == "transcript":
transcript = Transcript()
id = re.search(r'transcript_id "([^;]+)";?',attr).group(1)
parent = re.search(r'gene_id "([^;]+)";?', attr).group(1)
if not parent in list_gene:
continue
transcript.chr = chr
transcript.start = start
transcript.end = end
transcript.id = id
transcript.parent = parent
list_transcript[id] = transcript
elif type == "exon":
exon = Exon()
parent = re.search(r'transcript_id "([^;]+)";?',attr).group(1)
if not parent in list_transcript:
continue
exon.chr = chr
exon.start = start
exon.end = end
exon.parent = parent
list_exon.append(exon)
chr_gene(list_gene)
gene_len(list_gene)
gene_transcript(list_transcript)
transcript_exon(list_exon)
exon_pos(list_exon)
分别定义目标函数
chr_gene函数以记录了基因信息的字典作为输入，然后统计每条染色体上包含的基因数量。
def chr_gene(list_gene):
"""
param list_gene:
"""
print("Distribution of genes on genome")
count_gene = {}
for info in list_gene.values():
chr = info.chr
if chr in count_gene:
count_gene[info.chr] += 1
else:
count_gene[info.chr] = 1
with open("chr_gene.txt", 'w') as fp_out:
for chr,num in count_gene.items():
print("\t".join([chr,str(num)]) + "\n")
fp_out.write("\t".join([chr,str(num)]) +"\n")
gene_len函数以list_gene这个字典作为输入，统计每个基因的长度
def gene_len(list_gene):
"""
param list_gene:
"""
print("The length of gene")
with open("gene_len.txt",'w') as fp_out:
for gene_id,info in list_gene.items():
len = info.end - info.start +1
fp_out.write("\t".join([gene_id,str(len)]) + "\n")
print("\t".join([gene_id,str(len)]) + "\n")
gene_transcript函数以list_transcript作为输入，统计每个基因上的转录本的数量
def gene_transcript(list_transcript):
"""
:param list_transcript:
"""
print("The distributions of transcripts on genomes")
count_transcript = {}
for info in list_transcript.values():
gene_id = info.parent
if gene_id in count_transcript:
count_transcript[gene_id] += 1
else:
count_transcript[gene_id] =1
with open("gene_transcript.txt",'w') as fp_out:
for gene_id,num in count_transcript.items():
print("\t".join([gene_id,str(num)]) + "\n")
fp_out.write("\t".join([gene_id,str(num)]) + "\n")
transcript_exon 以list_exon作为输入，统计每个转录本上的外显子数量
def transcript_exon(list_exon):
"""
param list_exon
"""
print("Summarise of exons within transcript")
count_exon = {}
for exon in list_exon:
transcript_id = exon.parent
if transcript_id in count_exon:
count_exon[transcript_id] += 1
else:
count_exon[transcript_id] = 1
with open("transcript_exon.txt",'w') as fp_out:
for transcript_id, num in count_exon.items():
print("\t".join([transcript_id,str(num)]) + "\n")
fp_out.write("\t".join([transcript_id,str(num)]) + "\n")
exon_pos 以list_exon作为输入，统计转录本中每个外显子的位置信息
def exon_pos(list_exon):
"""
param list_exon
"""
print("The coordination of exons")
count_exon = {}
for exon in list_exon:
transcript_id = exon.parent
if transcript_id in count_exon:
count_exon[transcript_id] += ",%s-%s" % (str(exon.start),str(exon.end))
else:
count_exon[transcript_id] = "%s-%s" % (str(exon.start),str(exon.end))
with open("exon_pos.txt",'w') as fp_out:
for transcript_id,pos in count_exon.items():
print("\t".join([transcript_id,pos])+"\n")
fp_out.write("\t".join([transcript_id,pos]) + "\n")
调用主函数
if __name__ == "__main__":
main(args)