全面解析叶绿体基因组比较分析：必备软件与ggplot绘图脚本
前段时间，发表了一篇叶绿体基因组比较分析文章。尽管这类套路文已经不再热门，但在准备数据分析的过程中，仍感各种资料教程之零散，各种文章的做图工具也是令人眼花缭乱。故此，列出个人分析数据过程中的全套记录和工具选择，为后来者提供些微帮助。
1 叶绿体基因组的注释和可视化
注释使用曲小健老师的PGA，该软件可以批量注释多个组装好的叶绿体基因组。该软件有自带的注释很完整的参考基因组，可以满足大部分场合的基因组注释。不过也可以自己指定一个或者多个自己的参考基因组。
该软件可以在windows下本地安装运行，不过需要依赖Blastn+，简易但够用的参考命令如下：
PGA.pl -r plastomes-zhongzi -t target
-r 参考基因组的文件夹，可准备多个参考基因组，也可使用软件自带的参考基因组
-t 待注释基因组的文件夹，如果放入多个基因组fasta文件，即可批量注释
叶绿体基因组的可视化，一般使用OGDRAW在线工具，上传注释好的gb文件和配置好一些简单的参数，即可导出结果图，如下:
2 叶绿体基因组比较分析
2.1 长重复序列（repeat sequence analysis）分析
利用REPuter在线工具逐条基因组序列计算，大致的参考参数如下：
将多个序列的结果汇总成以下格式，下图是长重复序列的各类型个数：
下图是长重复的长度数据表格格式示意，每个样品各种重复的长度数据按列排布：
然后使用以下R语言脚本，即可一次成图，直接用于文章中：
library(ggplot2)
# 类型分布作图
# 读取数据
mydata_category <- read.table("category.txt", sep = "\t")
colnames(mydata_category) <- c("species","category","value") # 指定列名
mydata_category$category <- factor(mydata_category$category,levels = c("P","F","R","C")) # 指定序列类型的顺序，防止ggplot作图时自动排序
pic_category <- ggplot(mydata_category,mapping = aes(species,value,fill = category)) + # 分category填充颜色
geom_bar(stat='identity',position='dodge') +
# labs(x = 'species',y = 'time') +
# theme(axis.title =element_text(size = 16),axis.text =element_text(size = 14, color = 'black'))+
# theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
scale_x_discrete(labels = NULL) + # 不显示横座标
xlab(NULL) + ylab(NULL) + # 不显示横纵图示
theme(plot.margin = unit(c(0.3,0,0,1),"cm")) + # 指定画布边缘宽度，顺序为上，右，下，左
geom_text(label = mydata_category$value,vjust=-0.5,hjust=0.5, position = position_dodge(0.9),size = 2)
# 长度分布作图
mydata_longer <- read.table("longer.txt", sep = "\t",header = T,check.names = F)
# 统计长度数据分布
trans_mydata <- data.frame()
for(i in 1:length(colnames(mydata_longer))){
temp_data<- cbind(rep(colnames(mydata_longer)[i],4), c("30-40","40-50","50-60","60-70"),
table(cut(mydata_longer[,i],breaks = c(30,40,50,60,70),right = F)))
trans_mydata <- rbind(trans_mydata,temp_data)
}
trans_mydata[,3] <- as.numeric(as.character(trans_mydata[,3])) # 转换为数值
colnames(trans_mydata) <- c("species","length","value") # 指定列名
trans_mydata[,2] <- factor(trans_mydata[,2],levels = c("30-40","40-50","50-60","60-70")) #指定顺序，防止ggplot重排
pic_longer <- ggplot(trans_mydata,mapping = aes(species,value,fill = length)) + # 分length填充颜色
geom_bar(stat='identity',position='dodge') +
# labs(x = 'species',y = 'time') +
# theme(axis.title =element_text(size = 16),axis.text =element_text(size = 14, color = 'black'))+
theme(axis.text.x = element_text(angle = 30, hjust = 1)) + # 指定横坐标角度和对其位置
theme(plot.margin = unit(c(0,0,0,1),"cm")) + # 指定画布边缘宽度，顺序为上，右，下，左
xlab(NULL) + ylab(NULL) + # 不显示横纵图示
geom_text(label = trans_mydata$value,vjust=-0.5,hjust=0.5, position = position_dodge(0.9),size = 2)
# 两图合并
library(cowplot)
comb_pic <- plot_grid(pic_category,pic_longer,nrow = 2, align = "v", rel_heights = c(1,1.5),
labels = c("A","B"), label_x = 0, label_y = 1, hjust = -1, vjust = 1.5)
# 指定两图对齐方式和高度比例
# label_x = 1,label_y = 1,label_x = 0, label_y = 0, hjust = 0, vjust = 0 指定标签位置
comb_pic
ggsave("long-repeat-sequence.pdf",width = 10,height = 6) # 保存为pdf
成图的效果如下：
2.1 简单重复序列（SSR）分析
利用[Misa](Misa-web - IPK Gatersleben)在线工具计算SSR的重复次数、重复单元以及长度等各种数据，参数使用建议如下，当然也可参考其他发表文献：
使用fasta格式的为注释数据。该工具也有本地版本，为perl语言编写。
同长重复序列的分析类似，也是分为按类别（重复基序，motif）和重复长度来分类整理在线工具结果，如下：
然后使用脚本直接完成做图：
library(ggplot2)
# 长度分布作图
# 读取数据
mydata_times <- read.table("timesforr.txt", sep = "\t")
colnames(mydata_times) <- c("species","category","value") # 指定列名
mydata_times$category <- factor(mydata_times$category,levels = c("1","2","3","4","5")) # 指定序列类型的顺序，防止ggplot作图时自动排序
pic_times <- ggplot(mydata_times,mapping = aes(species,value,fill = category)) + # 分category填充颜色
geom_bar(stat='identity',position='dodge') +
# labs(x = 'species',y = 'time') +
# theme(axis.title =element_text(size = 16),axis.text =element_text(size = 14, color = 'black'))+
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
xlab(NULL) + ylab(NULL) + # 不显示横纵图示
scale_x_discrete(labels = NULL) + # 不显示横座标
geom_text(label = mydata_times$value,vjust=-0.5,hjust=0.5, position = position_dodge(0.9),size = 2) + # 柱型图上的数字
theme(legend.title = element_text(size = 8),legend.text = element_text(size = 8),
legend.key.size = unit(0.6,"lines")) + # 指定图例位置和字体大小
theme(plot.margin = unit(c(0.3,0,0,1),"cm")) # 指定画布边缘宽度，顺序为上，右，下，左
pic_times
# 类别分布作图
mydata_category <- read.table("category.txt", sep = "\t",check.names = F)
colnames(mydata_category) <- c("species","category","value") # 指定列名
mydata_category[,2] <- factor(mydata_category[,2],levels = c("A/T", "C/G","AG/CT", "AT/AT", "AAG/CTT", "AAT/ATT", "AAAG/CTTT",
"AAAT/ATTT", "AATT/AATT", "AGAT/ATCT","AAAAT/ATTTT", "AAAGG/CCTTT",
"AATTC/AATTG", "AAAAG/CTTTT")) #指定顺序，防止ggplot重排
pic_category <- ggplot(mydata_category,mapping = aes(species,value,fill = category)) + # 分longer填充颜色
geom_bar(stat='identity',position='dodge') +
# labs(x = 'species',y = 'time') +
# theme(axis.title =element_text(size = 16),axis.text =element_text(size = 14, color = 'black'))+
theme(axis.text.x = element_text(angle = 30, hjust = 1)) + # 指定横坐标角度和对齐位置
xlab(NULL) + ylab(NULL) + # 不显示横纵图示
geom_text(label = mydata_category$value,vjust=-0.5,hjust=0.5, position = position_dodge(0.9),size = 2.5) +
theme(legend.title = element_text(size = 8),legend.text = element_text(size = 8),legend.key.size = unit(0.6,"lines")) +
theme(plot.margin = unit(c(0,0,0,1),"cm")) # 指定画布边缘宽度，顺序为上，右，下，左
pic_category
# 两图合并
library(cowplot)
comb_pic <- plot_grid(pic_times,pic_category,nrow = 2, align = "v", rel_heights = c(1,1.5),
labels = c("A","B"), label_x = 0, label_y = 1, hjust = -1, vjust = 1.5) # 指定两图对齐方式和高度比例
# label_x = 1,label_y = 1,label_x = 0, label_y = 0, hjust = 0, vjust = 0 指定标签位置
comb_pic
ggsave("SSR.pdf",width = 12,height = 8) # 保存为pdf
出图效果如下：
2.3 密码子偏好
使用DAMBE逐个样品计算relative synonymous codon usage (RSCU)，然后汇总结果，由于本人在文章中直接提供表格，没有可视化呈现，此处就没有代码了。
2.4 IRs 边界变异分析
使用IRscope在线工具，上传gb注释文件，submit之后等待结果即可，但是服务器经常无响应，需要多次尝试。该工具输出结果仅为位图，无pdf或者svg等矢量图结果，且分辨率不高，宽度固定最高为2490像素，效果较为一般。
此外一次分析的序列上限为10个，更多数量的样品需要分多次完成。当然，多次结果的图片也需要自行Ps合并展示。之后，如发现更好工具，有必要淘汰IRscope。
论文中常见结果展示如下：
2.5 基因组变异展示
使用mVISIT在线工具，同时输入比对好的fasta文件和每个样品单独的注释文件，该注释文件格式特殊（如下图），可使用曲小健的perl脚本来生成。
然后对应上传到mVISTA进行运算：
提交后，结果会以网址的方式发送：
一般点选第一个VISTA-Point，或者选择任何一个，设置其为参考序列，不会出现在结果图中（下图）。因此，可以指定以下研究对象以外的类群作为参考基因组，一同比对后上传分析，该参考基因组不会出现在图中，或者也可以选择研究对象之一作为参考基因组，比对时就需要重复一遍该序列。
alignment PDF的结果如下，比对结果固定以四行显示：
但是该模式对序列较多时十分不友好，具体的图会被不随数量变化粗细的边框明显挡住：
解决方案有两个：
将多个序列用同一个参考序列分批计算，保证每次运行的序列数量较少，然后使用Adobe Illustrator等pdf编辑器工具拼接在一起，这是我现在选择的方式。每批次运算5-6个，可以较好平衡展示效果和拼图工作量。
-
不使用VISTA Point，直接使用PDF按钮导出结果：
结果如下，该方法图片效果更好，不会出现边框盖住内容的情况，但是固定以每行25K的长度展示，一个160k左右长度的序列会显示在7行（7页pdf）中，不利于文章展示。
此时仍然可以通过pdf编辑器拼接在一起，但是会出现最后一行较短的情况，略显遗憾。如https://doi.org/10.3389/fgene.2020.00802的做法，中间可以看到明显的拼接痕迹。
VISTA Browser尽管可以调整很多内容，但是当前基本处于不可用的状态，因为VISTA Browser要求的Java Applet插件在当前主流的浏览器中基本不可用了。
不过IE通过调整设置可以勉强运行，调整内容包括：
- Internet浏览器 > Internet选项 > 安全 > Internet > 自定义级别 > 把Active控件和脚本相关的选项打开
- 控制面板 > 程序 > Java（请保证已经安装好java）> 安全 > 编辑站点列表 > 将mVISTA服务器的网址https://pipeline.lbl.gov/添加进去，再次打开VISTA Browser之后经过一些确认即可进入。
但是进入后的页面仍然显示不完整，至此仍无法解决。
2.6 多态性位点分析
多态性位点需要逐个位点（包括基因和基因间隔区）计算其变异位点。具体步骤可能包括：
根据注释文件，将每个样品的每个位点序列单独提取出来
针对每个位点，合并多个样品的该位点序列
同时需要完成比对和变异位点计算
完成做图和可视化
这可能是一个相当繁琐的工作，所以多数文章使用了DNAsp滑移固定长度窗口的方式来展示计算并展示多态性位点，然后高变区域的位置倒退有哪些位点，如下图。但是，也应该注意到，此方式仍存在一些逻辑问题。另外这种取巧的方法就无法回避另一种常见的叶绿体基因组比较分析了，就是选择压力分析，因此此分析要求提取出所有cds的序列来计算Ka和Ks值了。
自此，本文提出一个个人的、半成品的但可行的流程，可以批量提取所有样品的各个位点，并批量完成比对和计算变异位点。
-
根据注释文件，将每个样品的每个位点序列单独提取出来
此步骤同样使用曲小健的perl脚本来完成叶绿体基因组中各位点。将所有注释好的gb文件放到同一个文件夹内，就可使用曲小健的脚本分别针对基因和基因间隔区进行批量提取。这些脚本的结果格式为同一个的gb文件的所有gene、cds或者基因间隔区会放到同一个fasta中。因此，多个fasta结果文件的相同位点还需要提取合并到一起才行。
-
针对每个位点，合并多个样品的该位点序列，并完成比对和变异参数计算
实际上由于叶绿体基因组的重复结构，不少位点都有两个拷贝，而且个别特别位点，比如rps12在分为三个区域且长度不一。因此，在合并不同提取结果中的同位点序列的时候还需要考虑拷贝问题。
因此，作者也花了很多精力写了一个R脚本（本文称为多态性批量计算脚本），可以针对曲小健脚本生成的多个提取结果，进行同基因合并、muscle多序列比对、变异参数计算（包括序列长度、比对gap数量、变异位点和变异比例、核苷酸多态性等）、基因名称规范（针对perl脚本的命名风格）的自动化实现，并生成一个csv汇总表格。本脚本还会同时输出每个位点的比对好的fasta文件，放置在同路径的fasta文件夹内，可用于后续可能的合并分析，比如针对类群关系较远的类群，仅用叶绿体基因组的cds区域或者genes来建树。脚本中默认不计算比对长度小于200bp的位点，也可以根据自己的需要对脚本进行修改。
具体使用上比较简单，将多个样品的perl脚本提取结果放到同一个文件夹中，然后在R中或者cmd中使用脚本即可，会实时显示运算结果。
# Date： 2022-08-05
# Author： Zhang Zhen
# Email： ficuszz@163.com
# Website: zhangzhen.plus
# Description： 读取当前文件夹的fasta文件，根据文件名提取出同源序列，批量比对，计算变异位点比例，然后写入scv文件
# fasta文件准备方式：使用曲小健博士的Perl脚本从叶绿体基因组gb文件中提取gene、CDS或者间隔区序列
# gb文件准备：下载叶绿体基因组全长序列，利用曲小健的PGA程序重新批量注释基因组，可保证所有注释结果风格类似，
# 各基因和片段名亦相同
# 载入所需库
library(stringr)
library(ape)
library(pegas)
# 第一部分，批量读取fasta文件内容
# 批量读取当前工作目录中的fasta文件，然后再根据fasta文件名整理出物种名列表，最后将文件中的fasta文件逐个复制给物种名列表
# 此处注意fasta文件的命名，必须统一为“genus_epithet_***.txt”
work_path <- choose.dir()
setwd(work_path)
path_list <- list.files(pattern = "fasta")
name_list <- list()
for(i in 1:length(path_list)){
name_list[[i]] <- read.FASTA(path_list[i])
}
# 在当前目录下新建fasta文件夹，用于存储后续的单基因比对文件
if (dir.exists("fasta")){
unlink("fasta", recursive = TRUE)
dir.create("fasta")
} else dir.create("fasta")
# 第二部分，指定一些序列名称和序列合并的函数
# 处理fasta文件中的序列名，删除perl脚本提取gene和intergenic时附带的物种名，
# 使得来自不同文件的同源基因名字相同，便于后续比较提取
seqname_deal <- function(x){
str_sub(x,1,str_locate("_", string = x)[1]-1)
}
# 相反的处理方式，去除序列Labels包含的序列名，仅含物种名，后续批量构建单基因树的时候可保证tip labels一致，便于比较
seqname_deal1 <- function(x){
str_sub(x,str_locate("_", string = x)[1]+1,nchar(x))
}
# 计算突变位点数量，变异位点和indel位点的综合
substi_number <- function(x){
ncol(x)-ncol(del.colgapsonly(x,threshold = 1/nrow(x)))+length(seg.sites(del.colgapsonly(x,threshold = 1/nrow(x))))
}
# 在第x个fasta文件中寻找第1个fasta文件中第i个基因的同名序列
combine_xto1 <- function(x){
for(j in 1:length(name_list[[x]])){
if(seqname_deal(as.alignment(name_list[[x]][j])$nam) == seqname_deal(as.alignment(name_list[[1]][i])$nam)){
temp_dna <<- c(temp_dna,name_list[[x]][j])
break;}
}
}
# 第三部分，处理IR区域的重复片段，且包括rps12反式剪切体
# 处理rps12的特殊情况，rps12在叶绿体上有三个区域,其中一个较短，约140bp，剩余两个对称位于IR区，此处删除掉较短的那个
# 避免对后续的多样性分析结果造成影响，此处筛选的标本是小于400bp，兼顾适当的长度变异
for(i in 1:length(name_list)){
# 为每个物种生成一个片段名列表
allseqname_lists <- unlist(lapply(labels(name_list[[i]]), seqname_deal))
# 如果某个rps12片段的长度短于400bp，则删除该片段
k_dellists <- c()
if("rps12" %in% allseqname_lists){
for(k in which("rps12" == allseqname_lists)){
if(length(name_list[[i]][[k]])<400){
k_dellists <- c(k_dellists,-k)
}
}
}
else {break;}
name_list[[i]] <- name_list[[i]][k_dellists]
}
# 处理重复序列，对所有的序列名称内的字符进行重排，重排后一致的为重复序列
reorder_seqname <- function(x){
paste(sort(unlist(strsplit(x,split = ""))), collapse = "")
}
allseqname_lists <- unlist(lapply(labels(name_list[[1]]), seqname_deal))
re_seqname <- unlist(lapply(allseqname_lists,reorder_seqname))
# 根据gene或者intergenic的不同类别 分别生成删除索引
delete.list <- c()
if(max(str_count(allseqname_lists, "-")) == 1){
for(i in 1:length(allseqname_lists)){
if(sum(allseqname_lists == allseqname_lists[i]) == 2 ){
delete.list <- c(delete.list, -which(allseqname_lists == allseqname_lists[i])[2])
}
}
} else {
for(i in 1:length(re_seqname)){
if(sum(re_seqname == re_seqname[i]) == 2 ){
delete.list <- c(delete.list, -which(re_seqname == re_seqname[i])[2])
}
}
}
# 执行删除
if(length(delete.list) != 0) {
name_list[[1]] <- name_list[[1]][unique(delete.list)]
}
# 第四部分，处理一些个人指定需要删除的序列
# 定义需要删除的结果，包括IR的反向重复基因，以及对含有内含子的基因的处理,如：
# trnK,rps16,trnG,ycf3,trnL,trnV,clpP,rpl16,trnI,trnA编码区过短，长度主要由内含子组成，在此不计入基因中
# ycf68为假基因，位于trnL的内含子中，在此仅计入trnL的内含子，不计算ycf68为基因
supple_gene <- c("trnK-UUU","rps16","trnG-UCC","ycf3","trnL-UAA","trnV-UAC",
"clpP","rpl16","trnI-GAU","trnA-UGC","ycf68")
allseqname_lists <- unlist(lapply(labels(name_list[[1]]), seqname_deal))
delete.lists <- -na.omit(match(supple_gene,allseqname_lists))
if(length(delete.lists) != 0){
name_list[[1]] <- name_list[[1]][delete.lists]
}
# 第五部分，程序主体，计算各类参数如比对长度，核苷酸多态性，gap数量
# 初始化
seqname_lists <- c()
seqlength_lists <- c()
rate_lists <- c()
sequence_number_list <- c()
varsite_lists <- c()
gaps_lists <- c()
segsite_lists <- c()
segsiterate_lists <- c()
# 执行程序主体，逐个计算同源序列的多样性参数
for(i in 1:length(name_list[[1]])){
temp_dna <- name_list[[1]][i]
# 在第2至最后一个fasta文件中，找第1个fasta文件第i个基因片算的同源序列，并写入变量中，待下一步计算
lapply(2:length(name_list),combine_xto1)
align_temp_dna <- muscle(temp_dna)
seqname_lists <- c(seqname_lists,seqname_deal(as.alignment(name_list[[1]][i])$nam))
# 将比对好的单基因fasta文件逐个写入fasta文件夹中
write.dna(updateLabel(align_temp_dna,labels(align_temp_dna),seqname_deal1(labels(align_temp_dna))),
paste("fasta/",as.character(seqname_deal(as.alignment(name_list[[1]][i])$nam)),".fas", sep = ""),format = "fasta")
# 逐个计算基因长度、变异位点数量、gap数量、总数量、比例和核苷酸多态性
seqlength_lists <- c(seqlength_lists,ncol(align_temp_dna))
varsite_lists <- c(varsite_lists, length(seg.sites(del.colgapsonly(align_temp_dna,threshold = 1/nrow(align_temp_dna)))))
gaps_lists <- c(gaps_lists, ncol(align_temp_dna)-ncol(del.colgapsonly(align_temp_dna,threshold = 1/nrow(align_temp_dna))))
sequence_number_list <- c(sequence_number_list, nrow(align_temp_dna))
segsite_lists <- c(segsite_lists,substi_number(align_temp_dna))
segsiterate_lists <- c(segsiterate_lists,substi_number(align_temp_dna)/ncol(align_temp_dna))
rate_lists <- c(rate_lists,nuc.div(align_temp_dna))
print(paste("第",i,"个片段",seqname_deal(as.alignment(name_list[[1]][i])$nam),
"的比对后序列长度为:",ncol(align_temp_dna),"核苷酸多样性为:",nuc.div(align_temp_dna),
"突变位点数量(含indel)为:",substi_number(align_temp_dna),
"突变位点数量比例为:",substi_number(align_temp_dna)/ncol(align_temp_dna)))
}
last_result <- t(rbind(seqname_lists,seqlength_lists,sequence_number_list,varsite_lists,gaps_lists,segsite_lists,segsiterate_lists,rate_lists))
# 第六部分，规则化输出结果，包括列名和基因名
colnames(last_result) <- c("lociname","alignment length","number of sequences","number of variable siets","number of gap","number of substitutionsite (including indels)",
"rates of substitution site (including indels)","nucleotide diversity Pi")
# 替换不规则序列名
replace_lists <- matrix(c("psbA-trnK-UUU-2","psbA-trnK-UUU","trnK-UUU-2-matK","trnK-UUU-matK","matK-trnK-UUU-1","matK-trnK-UUU",
"trnK-UUU-1-rps16-2","trnK-UUU-rps16","rps16-2-rps16-1","rps16-intron","rps16-1-trnQ-UUG","rps16-trnQ-UUG",
"trnS-GCU-trnG-UCC-1","trnS-GCU-trnG-UCC","trnG-UCC-1-trnG-UCC-2","trnG-UCC-intron",
"trnG-UCC-2-trnR-UCU","trnG-UCC-trnR-UCU","atpF-2-atpF-1","atpF-intron","atpF-1-atpH","atpF-atpH",
"rpoC1-2-rpoC1-1","rpoC1-intron","ycf3-3-ycf3-2","ycf3-intron-1","ycf3-2-ycf3-1","ycf3-intron-2",
"ycf3-1-trnS-GGA","ycf3-trnS-GGA","trnL-UAA-1-trnL-UAA-2","trnL-UAA-intron",
"trnL-UAA-2-trnF-GAA","trnL-UAA-trnF-GAA","ndhC-trnV-UAC-2","ndhC-trnV-UAC",
"trnV-UAC-2-trnV-UAC-1","trnV-UAC-intron","clpP-3-clpP-2","clpP-intron-1",
"clpP-2-clpP-1","clpP-intron-2","clpP-1-psbB","clpP-psbB","petB-1-petB-2","petB-intron",
"petB-2-petD-1","petB-petD","petD-1-petD-2","petD-intron","rpl16-2-rpl16-1","rpl16-intron",
"trnL-CAA-ndhB-2","trnL-CAA-ndhB","ndhB-1-rps7","ndhB-rps7","rps12-1-trnV-GAC","rps12-trnV-GAC",
"rrn16-trnI-GAU-1","rrn16-trnI-GAU","ndhA-2-ndhA-1","ndhA-intron","trnA-UGC-2-trnA-UGC-1","trnA-UGC-intron",
"trnI-GAU-1-rrn16","trnI-GAU-rrn16","trnK-UUU-1-rbcL","trnK-UUU-rbcL","trnV-UAC-1-trnV-UAC-2","trnV-UAC-intron",
"trnV-UAC-2-ndhC","trnV-UAC-ndhC","trnF-GAA-trnL-UAA-2","trnF-GAA-trnL-UAA",
"trnL-UAA-2-trnL-UAA-1","trnL-UAA-intron","trnL-UAA-1-trnT-UGU","trnL-UAA-trnT-UGU",
"trnS-GGA-ycf3-1","trnS-GGA-ycf3","ycf3-1-ycf3-2","ycf3-intron-2","ycf3-2-ycf3-3","ycf3-intron-1",
"ycf3-3-psaA","ycf3-psaA","rpoC1-1-rpoC1-2","rpoC1-intron","trnR-UCU-trnG-UCC-2","trnR-UCU-trnG-UCC",
"trnG-UCC-1-trnS-GCU","trnG-UCC-trnS-GCU","trnK-UUU-1-chlB","trnK-UUU-chlB",
"psaM-trnG-UCC-1","psaM-trnG-UCC","rpl2-2-rpl2-1","rpl2-intron","ndhA-1-ndhA-2","ndhA-intron",
"trnV-GAC-rps12-1","trnV-GAC-rps12","rps7-ndhB-1","rps7-ndhB","ndhB-1-ndhB-2","ndhB-intron",
"ndhB-2-trnL-CAA","ndhB-trnL-CAA","rpl2-1-rpl2-2","rpl2-intron","trnG-UCC-2-trnG-UCC-1","trnG-UCC-intron",
"atpH-atpF-1","atpH-atpF","atpF-1-atpF-2","atpF-intron","petD-2-rpoA","petD-rpoA",
"trnA-UGC-1-trnA-UGC-2","trnA-UGC-intron","trnG-UCC-1-trnT-GGU","trnG-UCC-trnT-GGU",
"trnI-GAU-2-trnI-GAU-1","trnI-GAU-intron","trnV-UAC-1-trnM-CAU","trnV-UAC-trnM-CAU",
"ndhB-2-ndhB-1","ndhB-intron","trnI-GAU-1-trnI-GAU-2","trnI-GAU-intron","rps12-2-rps12-1", "rps12-intron",
"trnT-UGU-trnL-UAA-1", "trnT-UGU-trnL-UAA","trnM-CAU-trnV-UAC-1", "trnM-CAU-trnV-UAC",
"rps12-1-rps12-2","rps12-1-rps12-2"),nrow = 2,ncol = 69)
for(i in 1:ncol(replace_lists)){
last_result[last_result == replace_lists[1,i]] <- replace_lists[2,i]
}
# 写入结果到文件
write.csv(last_result,"last_result.csv", quote = F, row.names = F)
print("请在当前工作目录下的last_result.csv文件中查看结果")
需声明，此脚本仍非常稚嫩，如有使用者，需要自行负责可靠性，可能也需要自行修改部分代码。
-
完成做图和可视化
输入结果稍作整理（格式如下），就可以使用以下脚本来完成。
脚本如下：
library(ggplot2)
library(cowplot)
library(showtext)
install.packages("showtext")
font_add('Arial','c:/Windows/Fonts/arial.ttf') # 解决字体加载错误问题
showtext_auto() # 自动调用showtext，否则无法在ggsave()中使用，因为ggsave会自动打开和关闭图形设备
gene <- read.table("gene.txt",header = T,sep = "\t")
# 按照出现顺选指定横座标顺序，预防下一步作图是ggplot自行排列顺序 https://www.jianshu.com/p/437906c6b063
gene$lociname <- factor(gene$lociname,levels = unique(gene$lociname))
# 全部数据作图
gene_plot <- ggplot(data = gene, aes(x = lociname, y = nucleotide.diversity.Pi, group = 1)) +
geom_bar(stat = 'identity',
fill = ifelse(gene$nucleotide.diversity.Pi > sort(gene$nucleotide.diversity.Pi,decreasing = T)[6],'#6bafd6','gray50')) +
# 条形图，并选择前5位位点标记独特颜色
geom_hline(yintercept = mean(gene$nucleotide.diversity.Pi)) + # 添加横线，横线位置为所有数值的平均值
annotate("text",x = "rrn16", y = 0.0019,label = "0.001676",size = 4, hjust = 0.65) + # 添加一个单独的数字到图中，概述中为横线的具体数值
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), text = element_text(size =14)) + #定义横座标字号和对其位置
labs(x=NULL,y=NULL) + # 不使用横纵坐标说明位子
geom_text(label = gene$alignment.length,vjust = 0.4, hjust = -0.3,
position = position_dodge(1.5),size = 4, angle = 90) + # 条形图上添加比对长度数值
scale_y_continuous(expand = expansion(mult = c(0.02,0.12))) # 调整灰色画布的上下边界，利于一些标注，使其不出框
gene_plot
help("scale_y_continuous")
intergeneric <- read.table("intergeneric.txt",header = T,sep = "\t",check.names = F)
intergeneric_plot <- ggplot(data = intergeneric, aes(x = lociname, y = `nucleotide diversity Pi`, group = 1)) +
geom_bar(stat = 'identity',
fill = ifelse(intergeneric$`nucleotide diversity Pi` > sort(intergeneric$`nucleotide diversity Pi`, decreasing = T)[6],'#6bafd6','gray50')) +
geom_hline(yintercept = mean(intergeneric$`nucleotide diversity Pi`)) + # 添加平均值水平线
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), text = element_text(10)) +
labs(x=NULL,y=NULL) +
geom_text(label = intergeneric$`alignment length`,vjust = 0.4, hjust = -0.3,
position = position_dodge(1.5),size = 3, angle = 90) + # 添加比对长度
annotate("text",x = "rrn4.5-rrn5", y = 0.0048,label = "0.004353",size = 4, hjust = 0.5) +
scale_y_continuous(expand = expansion(mult = c(0.02,0.08)))
intergeneric_plot
plot_grid(gene_plot, intergeneric_plot, nrow = 2, align = "v", rel_heights = c(0.75,1),
label_x = "name of loci", label_y = "nucleotide.diversity.Pi")
# 导出PDF文件
ggsave("test.pdf",width = 18,height = 10)
出图效果如下，包括基因和基因间隔区两种分类、比对长度、平均核苷酸多态性，以及多态性top5的蓝色突出显示。
2.7 选择压力分析
选择压力分析是作者之前基本没有涉及的领域，听说过的也就是CodeML、easycodeml以及datamonkey等，这些软件的理论较为复杂，而且也似乎不能批量的完成多位点多样品的选择压力分析，在叶绿体基因组比较分析中较难应用。而多个位点多样品的选择压力分析本身也较为复杂，需要以下步骤完成：
提取出注释好cds序列
合并多个样品的同一cds位点的序列并比对
针对任一个cds位点，计算参考序列和所有样品序列之间的Ka（非同义替换）和Ks（同义替换）两个参数，通过计算Ka/Ks来评估选择压力
批量完成所有cds位点的Ka/Ks的计算
可视化
对于第一步和第二步，和多态性位点分析中类似，同样先使用曲小健的perl脚本提取所有注释gb文件的cds序列，然后使用作者的多态性批量计算脚本进行逐个位点的cds序列合并并比对。
对于第三步和第四步的任务，需要考虑批量计算的可能行，作者选择了TBtools的simple Ka/Ks Calculator功能完成Ka/Ks的计算。该功能需要两个输入文件，一个是比对好的fasta文件，一个是序列名对的二列的表格，用来指定计算对象，因为每次Ka/Ks的计算都是针对两条序列的。注意，因为此步需要参考序列，所以在第一步和第二步那里，就需要将参考基因组一同运算在内。那么，批量计算的步骤就很明显了：
合并所有cds位点的fasta文件，也就是将作者脚本生成多个比对好的fasta格式的cds序列文件合并到一起。
-
生成一套Ka/Ks序列名称对，用以指导TBtools的计算。这套名称对，应该为每个cds位点，都指定好参考序列和所有的样品序列之间名称对。这个名称对的生成作者也利用R脚本生成，输入数据即为作者多态性批量计算脚本输出的fasta文件夹。该脚本也会同时合并所有的cds序列，可直接作为TBtools Ka/Ks计算的输入文件。注意，此脚本要求参考序列位于所有样品之后。
library(ape) # 批量获取文件夹内文件名 name_list <- list.files() # 批量读取文件夹内的全部ufasta文件至列表 fasta_list <- list() for(i in 1:length(list.files())){ fasta_list[[i]] <- read.FASTA(list.files()[i]) } name_pairs <- list() # 批量生成物种名称对 for(j in 1:length(fasta_list)){ # 给采集号为主的序列名加上位点名称后缀 names(fasta_list[[j]]) <- paste(names(fasta_list[[j]]),substr(name_list[j],1,nchar(name_list[1])-4),sep = "_") # 合并全部fasta文件 if(j==1){write.FASTA(fasta_list[[j]],"comb.txt")} else {write.FASTA(fasta_list[[j]],"comb.txt",append = TRUE)} # 按位点逐个生成序列名称对 name_for_one_gene <- names(fasta_list[[j]]) for(i in 1:(length(name_for_one_gene)-1)){ name_pairs[[i+12*j-12]] <- c(name_for_one_gene[i],tail(name_for_one_gene,1)) } # 此处的12为含参考序列的样品总数，需自行调整 } # 转置列表形式的序列名称对至矩阵形式 name_pairs_last <- t(matrix(unlist(name_pairs), ncol = length(name_pairs), nrow = 2)) # 写出至文件 write.table(name_pairs_last,"name_pairs_last.txt",quote = FALSE,sep = "\t",row.names = FALSE, col.names = FALSE) 将脚本生成序列对名称文件和合并好的fasta文件使用TBtools来计算Ka/Ks即可。
第五步可视化，而可视化涉及到一个问题。一般来说，Ka/Ks大于1，被认为是正向选择；小于1，反向选择或者纯化选择；等于1，中性选择。但是，在实际运算中，还包括一些异常值，这些异常值在近缘类群中也不算少见，但却在已发表文献中讨论较少。比如，Ka和Ks都是0的时候，软件会给出NaN值，这种一般也可以看作为中性选择，毕竟Ka和Ks都没有发生；还有就是Ka大于0而Ks等于0的时候，会出现无穷大的结果，尽管也可以认为是大于1的正向选择，但是在可视化中却难以处理。
此外，在叶绿体基因组比较分析中，常常由于多个样品，每个位点会有多个Ka/Ks值。而文献中通常仅展示一个值，可能为其平均值，但求其平均也会降低信息含量。又考虑到上段提到的异常值问题，作者给出一个较少见的箱形图可视化方案，同时展示多个Ka/Ks值的分布特征和异常值，具体如下。上图为无穷大异常值的数量分布，下图多个Ka/Ks值的箱形图，其中NaN值被视为1。
该图的生成需要准备三个数据文件和一个做图脚本。
数据文件1，剔除掉所有inf值的Ka/Ks结果表，NaN已被替换为1。该文件用来完成箱形图。
数据文件2，记录每个cds位点inf值出现次数的表格。该文件用来完成记录inf值出现次数的折线图。
数据文件3，是数据文件3的简化，仅记录inf值出现次数不为0的cds位点的结果。该文件用来标注折线图具体代表的出现次数，0次不标注。如果你对R语言做图非常熟悉，这个表格应为多余，可通过命令来完成0次不标注的效果。但作者R水平有限，采取另准备一个文件的笨办法。
而所需做图脚本如下：
library(ggplot2)
library(cowplot)
library(showtext)
font_add('Arial','c:/Windows/Fonts/arial.ttf') # 解决字体加载错误问题
showtext_auto()
# 读取文件
kaks_value <- read.table("kaks.txt",sep = "\t",header = T)
inf_fullvalue <- read.table("inf-allvalue.txt",sep = "\t",header = T)
inf_value <- read.table("inf.txt",sep = "\t",header = T)
# 作图
kaks <- ggplot(kaks_value,mapping = aes(gene, value)) +
geom_boxplot() +
labs(x=NULL,y=NULL) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), text = element_text(size = 15)) + # 调整横座标
# scale_y_continuous(expand = expansion(mult = c(0,0.5))) +
geom_line(data = inf_fullvalue, aes(gene,(value*0.05+1.6), group = 1)) +
annotate("text", x = inf_value$gene,y = (inf_value$value*0.05+1.65), label = inf_value$value) # 添加额外的备注（正无穷位点的值）
kaks
ggsave("kaks.pdf",width = 14,height = 5)
3. 后记和声明
本文记录所记录的分析工具和可视化方案，均为作者在发表文章中所用样式，仅为个人偏好，读者可酌情参考。
此外，尽管作者在本文中列举了诸多个人计算和绘图所用R脚本，但是其质量、效率、可靠性和通用性都不敢有所保证，仅为作者数据分析时所用工具。如需参考还需根据个人的工作路径、样品数量或者测序文件性质针对性的修改。所以可能仍需要读者具有一定的R语言基础，这些半成品脚本才有可参考的价值。因此，标题亦有夸张之嫌，仅为吸引读者，还请见谅！
上一篇： 理解并掌握 np.bincount 函数的使用方法
下一篇： 剖析质粒结构的学习过程 - 更新版