其他
m6A metagene distribution 计算详解
1引言
绘制 m6A peaks 在 mRNA
上的分布的工具有 metaplotR,Guitar 等,前者结合了 perl 语言和 R 语言来分析绘图,后者则是一个 R 包。
m6A 在 mRNA 上的分布在许多 m6A 相关的文章里也经常出现:
今天尝试全程使用 python 来进行分析和绘图。
大体思路:
筛选出 最长转录本 GTF 文件。 计算每个 mRNA 对应 5'UTR, CDS, 3'UTR 区域的长度。 将每个特征区域转换为 单个记录位点,与转录组位置对应起来。 将上步的数据与 m6A peaks 取交集然后 计算相对位置。
2获取最长转录本 GTF 文件
针对一个基因含有多个转录本,我们选取最长的转录本来分析:
import gzip
import pandas as pd
计算 mRNA 每个特征的长度
# 打开 gtf 文件
with gzip.open('hg19.ncbiRefSeq.gtf.gz','rt') as gtf:
# 信息保存在字典里
trans_len = {}
utr5_len = {}
cds_len = {}
for line in gtf:
# 跳过注释行
if line.startswith('#'):
continue
# 分割
fields = line.split()
# 类型
type = fields[2]
if type == 'exon':
# 名称
gene_name = fields[17].replace('"','').replace(';','')
gene_id = fields[9].replace('"','').replace(';','')
trans_id = fields[11].replace('"','').replace(';','')
# 连接名称
key = '|'.join([gene_name,gene_id,trans_id])
# 计算多个外显子长度
length = int(fields[4]) - int(fields[3]) + 1
# 累计求和
trans_len.setdefault(key,0)
trans_len[key] += length
elif type == 'CDS':
# 名称
gene_name = fields[17].replace('"','').replace(';','')
gene_id = fields[9].replace('"','').replace(';','')
trans_id = fields[11].replace('"','').replace(';','')
# 连接名称
key = '|'.join([gene_name,gene_id,trans_id])
# 计算多个CDS长度
length = int(fields[4]) - int(fields[3]) + 1
# 累计求和
cds_len.setdefault(key,0)
cds_len[key] += length
elif type == '5UTR':
# 名称
gene_name = fields[17].replace('"','').replace(';','')
gene_id = fields[9].replace('"','').replace(';','')
trans_id = fields[11].replace('"','').replace(';','')
# 连接名称
key = '|'.join([gene_name,gene_id,trans_id])
# 计算多个5'UTR长度
length = int(fields[4]) - int(fields[3]) + 1
# 累计求和
utr5_len.setdefault(key,0)
utr5_len[key] += length
else:
pass
else:
pass
整理为表格:
微信扫一扫付费阅读本文
可试读14%
微信扫一扫付费阅读本文