使用R语言计算微卫星重复序列单体的一致序列/序列两两差异碱基数

参考 TRASH 这个软件里的代码，利用TRASH软件的输出结果
运行的TRASH命令是
~/mingyan/biotools/TRASH-main/TRASH_run.sh ../at.col0.chr.fna –o trash.output –seqt at.Chr02.again/atCen178_templates.csv –par 16
输出文件夹和重复序列单体的模板文件路径用绝对路径
TRASH的代码里如何判断鉴定的重复序列单体属于模板里的单体这个还得仔细看看
随机挑选一些单体
library(tidyverse)
read_csv(“C:/Users/lenovo/Desktop/atGenomeScience/all.repeats.from.at.col0.chr.fna.csv”) %>%
filter(!is.na(class)) %>%
filter(width==178) %>%
sample_frac(0.01) %>%
mutate(id=row_number()) %>%
mutate(seqid=paste0(class,”_”,id)) -> temp.dat
将结果写出到fasta文件里
library(seqinr)
write.fasta(sequences = as.list(temp.dat$seq),
names = as.list(temp.dat$seqid),
file.out = “At.satellite.178.fa”,
open=”w”)
这个长度都是178bp，这里有一个疑问是 这个序列是否还需要去做比对，TRASH的代码里是做了比对，mafft的比对代码
mafft –quiet –inputorder –kimura 1 –retree 1 At.satellite.178.fa > At.satellite.178.aligned.fa
计算一致序列的代码
alignment<-read.alignment(“At.satellite.178.fa”,format = “fasta”)

consensus 这个函数来源于 seqinr 这个R包

consensus.sequence<-toupper(paste(consensus(alignment,method = “majority”),collapse = “”))
如果做多序列比对，会引入gap，计算一致性序列的话也会有gap,这个gap是否需要保留？
TRASH的代码里
consensus.sequence[consensus.sequence !=”-“]
我猜这个代码里是想把gap去掉，但是也起不到这个作用，暂时想不通这个代码是在做什么。
两条序列计算edit.distance
stringdist::stringdist(“ATCGAAAG”,”ATCGTTTT”,method = “lv”)
这里的方法是 Levenshtein distance 我查了一下这个方法如果两个序列的长度相等计算的就是两条序列差异的碱基数量
这个是计算一对序列的差异碱基数，如果要计算多条序列两两之间的差异可以用函数
Biostrings::stringDist(c(“lazy”, “HaZy”, “crAzY”),method = “levenshtein”)
介绍拟南芥着丝粒特征的论文
The genetic and epigenetic landscape of the Arabidopsis centromeres
https://www.science.org/doi/10.1126/science.abi7489
Fig2c 的热图数据应该可以用这个方法算

计算数据

read_csv(“C:/Users/lenovo/Desktop/atGenomeScience/all.repeats.from.at.col0.chr.fna.csv”) %>%
filter(seq.name==”Chr2″) %>%
filter(start>=4808990) %>%
filter(end<=4826790) %>%
pull(seq) %>%
Biostrings::stringDist(method = “levenshtein”) %>%
as.matrix() %>%
as.data.frame() %>%
rownames_to_column() %>%
pivot_longer(!rowname)
作图代码
heatmap.data %>%
magrittr::set_colnames(c(“x”,”y”,”value”)) %>%
mutate(x=factor(x,levels=1:100),
y=factor(y,levels=1:100)) %>%
#filter(x!=y) %>%
ggplot(aes(x=x,y=y))+
geom_tile(aes(fill=value))+
scale_fill_viridis_c()+
theme_bw(base_size = 15)+
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
legend.position = “top”,
legend.title = element_blank())+
labs(x=NULL,y=NULL)+
coord_equal()

有空专门出一期视频介绍作图代码
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！