附近含有 m6A 修饰的 Stop Codon 序列提取
附近含有 m6A 修饰的 Stop Codon 序列提取
背景
“m6A 修饰存在于 RNA 上，在 RNA 的出生，转运和成熟、降解、翻译等过程发挥着重要的作用。
”
m6A 修饰分布于 RNA 上的 5'UTR、CDS、3'UTR、Exon 和 Intron 里，主要富集在 stop codon 附近，那么 stop codon 附近的 m6A 修饰具体发挥着怎样的功能和分子机制研究的不是很清楚，在不同物种间的、不同组织间和不同生理状态下 m6A 修饰的水平也会随之动态变化。
AI 画个示意图如下：
m6A 位点的预测和实验研究验证的信息都有很多在线数据库可以查询，对于数据库里数据的分析可以是多样化和个性化的，讲述一个合理的生物学过程，结合实验验证，是文章的主题框架和思路。
人生和生命科学一样，在不断折腾和求真知的路途中找到真理和生命的本质！
准备
整合实验数据的 m6A 修饰位点注释信息的数据库有好几个，如 RMBase、RNAmod、REPIC 和 m6A-atals 等，可以查询单个基因的 m6A 位点信息或者直接下载全部的信息。此外还可查询下载其它 RNA 修饰的信息。
1、数据下载
直接去数据库下载，我下载了 m6A-atals 数据库的 Human 的 m6A 的 Basic_Site_Information 信息，下载方法是进入网址后点击 Download 进入下载界面，下载相应文件即可。
http://180.208.58.66/m6A-Atlas/index.html
或者 linux 里下载
$ wget http://180.208.58.66/m6A-Atlas/Download/m6A_Human_Basic_Site_Information.txt
下载好后的文件打开长这样：
包括染色体、起始位置、终止位置、分布区域、数据来源、实验类型等等
$ less -S m6A_Human_Basic_Site_Information.txt
chr start end pos strand name gene_id dist type seq data exper cell
human_m6A_1 chr1 564474 564474 1 + LOC101928626 ENSG00000225972 ncRNA_exonic unprocessed_pseudogene CAAGGTCAAAGGGAGTCCGAACTAGTCTCAGGCTTCAACAT No 25491922 GSE54921 PA-m6A-seq HeLa Control
human_m6A_2 chr1 564545 564545 1 + LOC101928626 ENSG00000225972 ncRNA_exonic unprocessed_pseudogene CTTCATAGCCGAATACACAAACATTATTATAATAAACACCC No 30867593 GSE121942 miCLIP HepG2 Control|25491922
human_m6A_3 chr1 564560 564560 1 + LOC101928626 ENSG00000225972 ncRNA_exonic unprocessed_pseudogene CACAAACATTATTATAATAAACACCCTCACCACTACAATCT No 30867593 GSE121942 miCLIP HepG2 SETD2 knockdown|25491922
然后在 excel 里打开对数据进行筛选，假如筛选细胞系 cell 为 CD8T、type 为 protein coding、dist 为 exonic，intronic 和 3UTR，然后提取 chr、start、end、strand、name、dist 列：
$ less m6a.txt
chr start end strand name dist
chr1 949587 949587 + ISG15 exonic
chr1 1168512 1168512 + B3GALT6 exonic
chr1 1169034 1169034 + B3GALT6 UTR3
2、提取 stop codon 信息
要提取 stop codon 信息首先得知道数据库用的注释文件或基因组版本是哪有个，不然序列位置和注释信息是对不上的。去 m6A-atals 找了一下也没看到有写，那就先下一个 hg38 的 gtf 文件看看，结果发现序列位置差的太远了，肯定不是，于是又下载了 hg19 的才对的上，毕竟早期的文献数据分析都是用的较早的基因组和注释文件，hg19 基因组是 2009 年出来的。
下载 UCSC hg19 注释文件：
$ wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/genes/hg38.ncbiRefSeq.gtf.gz
我们需要提取 stop codon 的位置信息：
chr、stop_codon、start、end、strand、gene_name 这么几列
$ zless -S hg19.ncbiRefSeq.gtf.gz | grep -w 'stop_codon' |sed 's/"/\t/g' |sed 's/;/\t/g' |awk '{print $1"\t"$3"\t"$4"\t"$5"\t"$7"\t"$18}' | wc -l
67419
有 67419 个 stop codon，一个基因有多个转录本，如果这多个转录本的 stop codon 位置不一样的话那么一个基因就会有多个不一样的 stop codon，但是这多个转录本的 stop codon 位置也可能是一样，所以这 67419 个终止密码子是有冗余或者重复的，我们需要去一下重复。
查看重复的 stop codon：
$ zless -S hg19.ncbiRefSeq.gtf.gz | grep -w 'stop_codon' |sed 's/"/\t/g' |sed 's/;/\t/g' |awk '{print $1"\t"$3"\t"$4"\t"$5"\t"$7"\t"$18}' |less
chr22 stop_codon 38689293 38689295 - CSNK1E
chr22 stop_codon 38689293 38689295 - CSNK1E
chr22 stop_codon 38689293 38689295 - TPTEP2-CSNK1E
chr22 stop_codon 38617476 38617478 - TMEM184B
chr22 stop_codon 38617476 38617478 - TMEM184B
chr22 stop_codon 38617476 38617478 - TMEM184B
chr22 stop_codon 38610883 38610885 + MAFF
chr22 stop_codon 38610883 38610885 + MAFF
chr22 stop_codon 38610883 38610885 + MAFF
chr22 stop_codon 38610883 38610885 + MAFF
去除重复：
$ zless -S hg19.ncbiRefSeq.gtf.gz | grep -w 'stop_codon' |sed 's/"/\t/g' |sed 's/;/\t/g' |awk '{print $1"\t"$3"\t"$4"\t"$5"\t"$7"\t"$18}' |sort -k3,3n |uniq |wc -l
27734
# 输出到文件保存结果
$ zless -S hg19.ncbiRefSeq.gtf.gz | grep -w 'stop_codon' |sed 's/"/\t/g' |sed 's/;/\t/g' |awk '{print $1"\t"$3"\t"$4"\t"$5"\t"$7"\t"$18}' |sort -k3,3n |uniq > stopc.txt
最后我们把 stopc.txt 第二列放到最后一列，给它加上和 m6a.txt 一样的列名就可以了。
$ less stopc.txt
chr start end strand name dist
chr1 70006 70008 + OR4F5 stop_codon
chr1 368595 368597 + OR4F29 stop_codon
chr1 879531 879533 + SAMD11 stop_codon
算法思路
计算每个 m6A 修饰位点距离 stop codon 的距离，满足一定距离（100bp）的保存对应的 stop codon 位置信息。
分析：
1、m6A 修饰位点可位于 stop codon 上游也能是下游
2、一个基因存在多个 stop codon 位点 （一样位置的stop codon已经去除）
3、一个基因存在多个 m6A 修饰位点
YTHDF3 有 8 个转录本：
要是一个个去 excel 里找到每个基因对应的 m6A 位点和 stop codon 位点再互相减，挑选出在 100bp 以内的 stop codon 位点信息。如果有多个 stop codon 得更加多次的相减，几万个基因得减到猴年马月！此时对于这样一次次条件判断的重复过程我们可以使用编程的循环语句就能轻松实现。
解决思路
1、循环每个基因，取出相应的 m6A 和 stop codon 信息
2、判断该基因的 stop codon 的数量
3、如果 stop codon ==1，则与每个 m6A 位置相减，取出绝对值小于 100bp 的 stop codon
3、如果 stop codon >1，循环该基因每个 stop codon，与每个 m6A 位置相减，取出绝对值小于 100bp 的 stop codon
和师妹共同讨论了文件和解决思路，一起完成了一个循环。以下是代码部分：
# 加载R包
library(tidyverse)
# 1、读取筛选过后的m6a文件和stop codon信息文件
stopcodon <- read.delim("C:\\Users\\admin\\Desktop\\stopc.txt")
m6A_293 <- read.delim("C:\\Users\\admin\\Desktop\\m6a.txt")
# 2、取共有基因交集
m6A_gene <- unique(m6A_293$name)
stop_gene <- unique(stopcodon$name)
oveloap <- intersect(m6A_gene,stop_gene)
# 3、去除没有m6a或者stop codon信息的基因
m6A_gene_clean <- m6A_293[which(m6A_293$name %in% oveloap ),]
stopcodon_clean <-stopcodon[which(stopcodon$name %in% oveloap),]
# 4、在stop codon clean信息文件添加一列，为终止密码子中间位置
stopcodon_clean$pos <- stopcodon_clean$start+1
# 5、循环筛选符合m6a位点在密码子上下游15nt 以内符合条件的stop codon信息
for(i in oveloap){
temp_m6A <- m6A_gene_clean %>% filter(name == i)
temp_stop <- stopcodon_clean %>% filter(name == i)
if(length(temp_stop$pos) == 1){
test_16 <- abs(temp_m6A$start - temp_stop$pos) <= 100
if("TRUE" %in% test_16){
tar_stop1 <- temp_stop
file1 <- paste(i,'s_tar_stop',sep="_")
assign(file1,temp_stop)
}else {
}
}else{
for(m in temp_stop$pos){
test_16m <- abs(temp_m6A$start - m) <= 100
# print(test_16m)
# print(i)
if("TRUE" %in% test_16m){
tar_stop2 <- temp_stop %>% filter(pos == m)
file2 <- paste(i,'m_tar_stop',sep="_")
assign(file2,tar_stop2)
}
else{
}
}
}
}
# 6、查看含有单个密码子和多个密码子的筛选结果数量
length(ls(pattern = '*s_tar_stop'))
[1] 725
length(ls(pattern = '*m_tar_stop'))
[1] 279
# 7、储存总的结果变量名
res_name <- ls(pattern ='*_[s,m]_tar_stop' )
# 8、结果数量
length(res_name)
[1] 1004
# 9、合并所有结果到一个数据表里
nm <- map(res_name,as.name)
final_res <- do.call('rbind',nm)
# 10、整理成bed格式的文件便于提取序列
final_res_bed <- final_res[,c(1,2,3,5,6,4,7)]
# 11、提取stop codon序列的文件
final_res_bed2fa <- final_res_bed
final_res_bed2fa$start <- final_res_bed2fa$start-1
head(final_res_bed)
chr start end name dist strand pos
1 chr10 101611386 101611388 ABCC2 stop_codon + 101611387
2 chr21 43716500 43716502 ABCG1 stop_codon + 43716501
3 chr3 52015032 52015034 ABHD14A stop_codon + 52015033
4 chr6 24701841 24701843 ACOT13 stop_codon + 24701842
5 chr4 2930248 2930250 ADD1 stop_codon + 2930249
6 chr2 228419209 228419211 AGFG1 stop_codon + 228419210
# 12、保存结果文件
write.table(final_res,file="C:\\Users\\admin\\Desktop\\final_res.csv",quote = F,sep=",",row.names = F)
write.table(final_res_bed2fa,file="C:\\Users\\admin\\Desktop\\final_res_bed2fa.csv",quote = F,sep=",",row.names = F)
这里产生的结果文件会以单个文件保存在环境变量中，环境变量会生成很多文件名。
代码优化：
想了一下其实不用分 stop codon 的数量，不管有几个，直接与每个 m6A 位置相减，取出绝对值小于 100bp 的 stop codon 就行了，其次把每次循环产生的结果保存追加到一个 list 中，这样环境变量就不会生成很多文件了。
以下是优化后的代码部分，循环前的代码是一样的，就不放了，这里放循环及后面代码：
可试读54%