m6A enrichment on peak center
测试开头
点击上方，关注老俊俊！
测试结尾
1引言
Transcriptome-wide mapping of N6-methyladenosine by m6A-seq based on immunocapturing and massively parallel sequencing
此文章包含具体数据分析代码。
如下图:
计算具体 某个 motif 里 peak 中点的相对距离的图,上图为 m6A 基序里 peak summit 的相对距离。m6A 基序通常为 RRACH(R(A/G),H(A/CU)) 。
当然这不仅局限于 m6A 的基序,还可以拓展到其它的 motif 等等。
思路比较简单,即是计算每个 motif 离 peak 中点的相对距离。首先我们需要获取 peak 的序列,然后确定我们需要找的 motif 即可。 稍微麻烦的是具体细节的处理。
接下来我已经写成了函数,来解读分享一下。
2提取序列
导入相关库:
import re
import seaborn as sns
import matplotlib.pyplot as plt
from pyfaidx import Fasta
提取 peaks 序列函数
def GetPeaksFa(peak,genomefie,outfasta,type='bed',a=0,b=0,minlen=5):
# load geneome
genome = Fasta(genomefie)
outfa = open(outfasta,'w')
# process peaks file
with open(peak,'r') as bed:
for line in bed:
fields = line.replace('n','').split()
chr = fields[0]
strand = fields[5]
start = int(fields[1])
end = int(fields[2])
if type == 'bed':
end = end
elif type == 'narrowpeak':
end = end + 1
else:
print('please mark your peaks file type(bed/narrowpeak)')
break
# extend upstream and downstram peaks
if strand == '+':
seq = genome[chr][(start - a):(end + b)].seq
elif strand == '-':
seq = genome[chr][(start + a):(end -b)].complement.reverse.seq
else:
seq = genome[chr][(start - a):(end + b)].seq
# mimimum seq length
if len(seq) >= minlen:
# fa name
name = fields[3] + '|' + '|'.join([chr,str(start),str(end)])
# output seq
outfa.write('>' + name + '|' + strand + 'n')
outfa.write(seq + 'n')
outfa.close()
如果你的 bed 文件
第六列包含 正确的正负链信息(+/-),则按照正负链进行序列提取,若不是则按照正链的规则提取。
参数解释:
peak: peaks 文件。 genomefie: 基因组序列文件。 outfasta: 输出序列文件名。 type: peaks 文件类型,分为 bed和narrowpeak两种类型,默认是 bed,因为这两种文件坐标类型不一样,前者为 0-based,后者为 1-based,所以计算会有差异。a/b: 对 peak 文件上下游进行拓展多少碱基序列,例如使用的单个碱基的 bed 文件如 macs2 的 summit结果,则可以使用该参数进行拓展,默认为 0。minlen: 指定最短序列的长度,小于这个阈值则丢弃,我这里用 m6A(RRACH motif)则默认指定了 5 个碱基最短。
使用提取序列:
GetPeaksFa(peak='GSE110320_HepG2_shCtrl_m6a_rep1-3_peak.bed',
genomefie='hg19.fa',type='bed',
outfasta='m6A-peaks_sequnce.fa')
以上是提取到的序列。
3计算 motif 离 peak center 的距离
定义函数
def CalculatePeaksRelativeDist(peaksfa,motif='[G|A][G|A]AC[A|C|T]|[T|G|A]GT[C|T][C|T]',mid=3):
# save in list
relpos = []
# open fa
with open(peaksfa,'r') as seq:
for line in seq:
line = line.replace('n','')
if line.startswith('>'):
next
peakmid = round(len(line)/2,ndigits=0)
# search all motif
pattern = re.compile(motif,flags=re.IGNORECASE)
pos_res = pattern.finditer(line)
# shift to A site position
allpos = [pos.start() + mid for pos in pos_res]
# calculate relative to peak center bases
relposition = [pos - peakmid for pos in allpos]
# save in list
relpos.extend(relposition)
return relpos
相对距离若在中点前面则为负数,若在后面则为正数。
参数解释:
peaksfa: peaks 对应的 fasta 序列文件。 motif: 你的 motif 序列,支持正则表达式, 忽略大小写,默认是所有 RRACH 的序列组合,包括反向序列,这里包括了会匹配负链的却按照正链提取序列的情况, 正则表达式为 [G|A][G|A]AC[A|C|T]|[T|G|A]GT[C|T][C|T],如果你的 bed 第六列是正确的正负链信息,你也可以使用[G|A][G|A]AC[A|C|T]来完全匹配,排除前面那种情况,你还可以直接指定具体的 motif 如GGACT等。mid: 指定计算的 motif 中点, 默认为 3,即为 RR ACH 的 A 为计算起始。
计算,这里计算了三个类型的 m6A:
GGACT = CalculatePeaksRelativeDist(peaksfa='m6A-peaks_sequnce.fa',
motif='GGACT',mid=3)
AAACT = CalculatePeaksRelativeDist(peaksfa='m6A-peaks_sequnce.fa',
motif='AAACT',mid=3)
GAACC = CalculatePeaksRelativeDist(peaksfa='m6A-peaks_sequnce.fa',
motif='GAACC',mid=3)
返回的结果为三个列表。
4可视化
plt.figure(figsize=[22,6])
sns.set(style='ticks')
plt.subplot(1,2,1)
# label size
plt.tick_params(labelsize=18)
# density plot
sns.kdeplot(data=GGACT,shade=False,color='red',alpha=1,legend=True,bw_adjust=0.005)
sns.kdeplot(data=AAACT,shade=False,color='green',alpha=1,legend=True,bw_adjust=0.005)
sns.kdeplot(data=GAACC,shade=False,color='blue',alpha=1,legend=True,bw_adjust=0.005)
# labels
plt.xlabel("m6A(RRACH) motif relative to peak center",fontsize=20)
plt.ylabel("peak density",fontsize=20)
# x lines
plt.axvline(x=0,ls="--",lw=2,c='black')
plt.xlim(-10000,10000)
plt.legend(['GGACT','AAACT','GAACC'],frameon=False,fontsize = 16,loc='upper left')
plt.subplot(1,2,2)
# label size
plt.tick_params(labelsize=18)
# density plot
sns.kdeplot(data=[i for i in GGACT if i >= -3000 and i <= 3000],shade=False,color='#DA1212',alpha=1,label = 'test',bw_adjust=0.1)
# labels
plt.xlabel("m6A(RRACH) motif relative to peak center",fontsize=20)
plt.ylabel("peak density",fontsize=20)
# x lines
plt.axvline(x=0,ls="--",lw=2,c='black')
plt.xlim(-3000,3000)
plt.legend(['GGACT'],frameon=False,fontsize = 16,loc='upper left')
5结果验证
下面是一个参考结果:
重新计算绘图:
relpos1 = CalculatePeaksRelativeDist(peaksfa='Bouleall.tag.uniq.clean.CITS.001.singletion.3UTRTTS.101nt.fa',
motif='GTTT',mid=0)
relpos2 = CalculatePeaksRelativeDist(peaksfa='Bouleall.tag.uniq.clean.CITS.001.singletion.3UTRTTS.101nt.fa',
motif='TTT',mid=0)
relpos3 = CalculatePeaksRelativeDist(peaksfa='Bouleall.tag.uniq.clean.CITS.001.singletion.3UTRTTS.101nt.fa',
motif='GTT',mid=0)
# plot
plt.figure(figsize=[8,5])
sns.set(style='ticks')
# plt.subplot(1,2,1)
# label size
plt.tick_params(labelsize=18)
# density plot
sns.kdeplot(data=relpos1,shade=False,color='#DA1212',alpha=1,bw_adjust=0.09)
sns.kdeplot(data=relpos2,shade=False,color='purple',alpha=1,bw_adjust=0.09)
sns.kdeplot(data=relpos3,shade=False,color='green',alpha=1,bw_adjust=0.09)
# labels
plt.xlabel("(GUUU) motif relative to peak center",fontsize=20)
plt.ylabel("peak density",fontsize=20)
# x lines
plt.axvline(x=0,ls="--",lw=2,c='black')
plt.axvline(x=-2,ls="--",lw=2,c='grey')
plt.xlim(-20,20)
plt.xticks(ticks=[-20,-10,-6,-7,-8,-9,-5,-4,-3,-2,-1,0,10,20],labels=[-20,-10,'','','','','','','','','',0,10,20])
plt.legend(['GUUU','UUU','GUU'],frameon=False,fontsize = 16,loc='upper left')
6结尾
最后面这部分数据不提供分享。
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦,数据代码已上传至QQ群,欢迎加入下载。
群二维码：
老俊俊微信：
知识星球：
所以今天你学习了吗？
今天的分享就到这里了，敬请期待下一篇！
最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！
如果觉得对您帮助很大，赏杯快乐水喝喝吧！
往期回顾
◀m6A metagene distribution 纠正坐标
◀m6A metagene distribution 计算详解
◀跟着 nature cell biology 学绘图–小提琴图
◀...
请关注“恒诺新知”微信公众号，感谢“R语言“，”数据那些事儿“，”老俊俊的生信笔记“，”冷🈚️思“，“珞珈R”，“生信星球”的支持！
