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