前段時間贯涎，發(fā)表了一篇葉綠體基因組比較分析文章伯铣。盡管這類套路文已經(jīng)不再熱門斗这，但在準(zhǔn)備數(shù)據(jù)分析的過程中，仍感各種資料教程之零散捎稚，各種文章的做圖工具也是令人眼花繚亂。故此求橄，列出個人分析數(shù)據(jù)過程中的全套記錄和工具選擇今野，為后來者提供些微幫助。

1 葉綠體基因組的注釋和可視化

注釋使用曲小健老師的PGA罐农，該軟件可以批量注釋多個組裝好的葉綠體基因組条霜。該軟件有自帶的注釋很完整的參考基因組，可以滿足大部分場合的基因組注釋啃匿。不過也可以自己指定一個或者多個自己的參考基因組蛔外。

該軟件可以在windows下本地安裝運(yùn)行，不過需要依賴Blastn+溯乒，簡易但夠用的參考命令如下：

PGA.pl -r plastomes-zhongzi -t target

-r 參考基因組的文件夾夹厌，可準(zhǔn)備多個參考基因組，也可使用軟件自帶的參考基因組
-t 待注釋基因組的文件夾裆悄，如果放入多個基因組fasta文件矛纹，即可批量注釋

葉綠體基因組的可視化，一般使用OGDRAW在線工具光稼，上傳注釋好的gb文件和配置好一些簡單的參數(shù)或南，即可導(dǎo)出結(jié)果圖，如下:

2 葉綠體基因組比較分析

2.1 長重復(fù)序列（repeat sequence analysis）分析

利用REPuter在線工具逐條基因組序列計(jì)算艾君，大致的參考參數(shù)如下：

將多個序列的結(jié)果匯總成以下格式采够，下圖是長重復(fù)序列的各類型個數(shù)：

下圖是長重復(fù)的長度數(shù)據(jù)表格格式示意，每個樣品各種重復(fù)的長度數(shù)據(jù)按列排布：

然后使用以下R語言腳本冰垄，即可一次成圖蹬癌，直接用于文章中：

library(ggplot2)
# 類型分布作圖
# 讀取數(shù)據(jù)
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) + # 不顯示橫座標(biāo)
  xlab(NULL) + ylab(NULL) + # 不顯示橫縱圖示 
  theme(plot.margin = unit(c(0.3,0,0,1),"cm")) + # 指定畫布邊緣寬度虹茶，順序?yàn)樯鲜判剑遥潞铮?  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)
# 統(tǒng)計(jì)長度數(shù)據(jù)分布
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]))  # 轉(zhuǎn)換為數(shù)值
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)) + # 指定橫坐標(biāo)角度和對其位置
  theme(plot.margin = unit(c(0,0,0,1),"cm")) + # 指定畫布邊緣寬度，順序?yàn)樯弦牛衣采觯吕。?  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 指定標(biāo)簽位置
comb_pic
ggsave("long-repeat-sequence.pdf",width = 10,height = 6) # 保存為pdf

成圖的效果如下：

2.1 簡單重復(fù)序列（SSR）分析

利用[Misa](Misa-web - IPK Gatersleben)在線工具計(jì)算SSR的重復(fù)次數(shù)、重復(fù)單元以及長度等各種數(shù)據(jù)封豪，參數(shù)使用建議如下崖瞭，當(dāng)然也可參考其他發(fā)表文獻(xiàn)：

使用fasta格式的為注釋數(shù)據(jù)。該工具也有本地版本撑毛，為perl語言編寫书聚。

同長重復(fù)序列的分析類似，也是分為按類別（重復(fù)基序藻雌，motif）和重復(fù)長度來分類整理在線工具結(jié)果雌续，如下：

然后使用腳本直接完成做圖：

library(ggplot2)

# 長度分布作圖
# 讀取數(shù)據(jù)
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) + # 不顯示橫座標(biāo)
  geom_text(label = mydata_times$value,vjust=-0.5,hjust=0.5, position = position_dodge(0.9),size = 2) + # 柱型圖上的數(shù)字
  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")) # 指定畫布邊緣寬度胯杭，順序?yàn)樯涎倍牛遥伦龈觯?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)) + # 指定橫坐標(biāo)角度和對齊位置
  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")) # 指定畫布邊緣寬度，順序?yàn)樯暇优彝缙担拢?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 指定標(biāo)簽位置
comb_pic
ggsave("SSR.pdf",width = 12,height = 8) # 保存為pdf

出圖效果如下：

2.3 密碼子偏好

使用DAMBE逐個樣品計(jì)算relative synonymous codon usage (RSCU)太闺，然后匯總結(jié)果糯景，由于本人在文章中直接提供表格，沒有可視化呈現(xiàn)省骂，此處就沒有代碼了蟀淮。

2.4 IRs 邊界變異分析

使用IRscope在線工具，上傳gb注釋文件钞澳，submit之后等待結(jié)果即可怠惶，但是服務(wù)器經(jīng)常無響應(yīng)，需要多次嘗試轧粟。該工具輸出結(jié)果僅為位圖策治，無pdf或者svg等矢量圖結(jié)果，且分辨率不高逃延，寬度固定最高為2490像素览妖，效果較為一般轧拄。

此外一次分析的序列上限為10個揽祥，更多數(shù)量的樣品需要分多次完成。當(dāng)然檩电，多次結(jié)果的圖片也需要自行Ps合并展示拄丰。之后府树，如發(fā)現(xiàn)更好工具，有必要淘汰IRscope料按。

論文中常見結(jié)果展示如下：

2.5 基因組變異展示

使用mVISIT在線工具奄侠，同時輸入比對好的fasta文件和每個樣品單獨(dú)的注釋文件，該注釋文件格式特殊（如下圖）载矿，可使用曲小健的perl腳本來生成垄潮。

然后對應(yīng)上傳到mVISTA進(jìn)行運(yùn)算：

提交后，結(jié)果會以網(wǎng)址的方式發(fā)送：

一般點(diǎn)選第一個VISTA-Point闷盔，或者選擇任何一個弯洗，設(shè)置其為參考序列，不會出現(xiàn)在結(jié)果圖中（下圖）逢勾。因此牡整，可以指定以下研究對象以外的類群作為參考基因組，一同比對后上傳分析溺拱，該參考基因組不會出現(xiàn)在圖中逃贝，或者也可以選擇研究對象之一作為參考基因組，比對時就需要重復(fù)一遍該序列迫摔。

alignment PDF的結(jié)果如下沐扳，比對結(jié)果固定以四行顯示：

但是該模式對序列較多時十分不友好，具體的圖會被不隨數(shù)量變化粗細(xì)的邊框明顯擋拙湔肌：

解決方案有兩個：

1. 將多個序列用同一個參考序列分批計(jì)算迫皱，保證每次運(yùn)行的序列數(shù)量較少，然后使用Adobe Illustrator等pdf編輯器工具拼接在一起辖众，這是我現(xiàn)在選擇的方式卓起。每批次運(yùn)算5-6個，可以較好平衡展示效果和拼圖工作量凹炸。

2. 不使用VISTA Point戏阅，直接使用PDF按鈕導(dǎo)出結(jié)果：

   
   

結(jié)果如下，該方法圖片效果更好啤它，不會出現(xiàn)邊框蓋住內(nèi)容的情況奕筐，但是固定以每行25K的長度展示，一個160k左右長度的序列會顯示在7行（7頁pdf）中变骡，不利于文章展示离赫。

此時仍然可以通過pdf編輯器拼接在一起，但是會出現(xiàn)最后一行較短的情況塌碌，略顯遺憾渊胸。如https://doi.org/10.3389/fgene.2020.00802的做法，中間可以看到明顯的拼接痕跡台妆。

VISTA Browser盡管可以調(diào)整很多內(nèi)容翎猛，但是當(dāng)前基本處于不可用的狀態(tài)胖翰，因?yàn)閂ISTA Browser要求的Java Applet插件在當(dāng)前主流的瀏覽器中基本不可用了。

不過IE通過調(diào)整設(shè)置可以勉強(qiáng)運(yùn)行切厘，調(diào)整內(nèi)容包括：

1. Internet瀏覽器 > Internet選項(xiàng) > 安全 > Internet > 自定義級別 > 把Active控件和腳本相關(guān)的選項(xiàng)打開

2. 控制面板 > 程序 > Java（請保證已經(jīng)安裝好java）> 安全 > 編輯站點(diǎn)列表 > 將mVISTA服務(wù)器的網(wǎng)址https://pipeline.lbl.gov/添加進(jìn)去萨咳，再次打開VISTA Browser之后經(jīng)過一些確認(rèn)即可進(jìn)入。

但是進(jìn)入后的頁面仍然顯示不完整疫稿，至此仍無法解決培他。

2.6 多態(tài)性位點(diǎn)分析

多態(tài)性位點(diǎn)需要逐個位點(diǎn)（包括基因和基因間隔區(qū)）計(jì)算其變異位點(diǎn)。具體步驟可能包括：

1. 根據(jù)注釋文件遗座，將每個樣品的每個位點(diǎn)序列單獨(dú)提取出來

2. 針對每個位點(diǎn)靶壮，合并多個樣品的該位點(diǎn)序列

3. 同時需要完成比對和變異位點(diǎn)計(jì)算

4. 完成做圖和可視化

這可能是一個相當(dāng)繁瑣的工作，所以多數(shù)文章使用了DNAsp滑移固定長度窗口的方式來展示計(jì)算并展示多態(tài)性位點(diǎn)员萍，然后高變區(qū)域的位置倒退有哪些位點(diǎn)腾降，如下圖。但是碎绎，也應(yīng)該注意到螃壤，此方式仍存在一些邏輯問題。另外這種取巧的方法就無法回避另一種常見的葉綠體基因組比較分析了筋帖，就是選擇壓力分析奸晴，因此此分析要求提取出所有cds的序列來計(jì)算Ka和Ks值了。

自此日麸，本文提出一個個人的寄啼、半成品的但可行的流程，可以批量提取所有樣品的各個位點(diǎn)代箭，并批量完成比對和計(jì)算變異位點(diǎn)墩划。

1. 根據(jù)注釋文件，將每個樣品的每個位點(diǎn)序列單獨(dú)提取出來

   此步驟同樣使用曲小健的perl腳本來完成葉綠體基因組中各位點(diǎn)嗡综。將所有注釋好的gb文件放到同一個文件夾內(nèi)乙帮，就可使用曲小健的腳本分別針對基因和基因間隔區(qū)進(jìn)行批量提取。這些腳本的結(jié)果格式為同一個的gb文件的所有g(shù)ene极景、cds或者基因間隔區(qū)會放到同一個fasta中察净。因此，多個fasta結(jié)果文件的相同位點(diǎn)還需要提取合并到一起才行盼樟。

2. 針對每個位點(diǎn)氢卡，合并多個樣品的該位點(diǎn)序列，并完成比對和變異參數(shù)計(jì)算

   實(shí)際上由于葉綠體基因組的重復(fù)結(jié)構(gòu)晨缴，不少位點(diǎn)都有兩個拷貝译秦，而且個別特別位點(diǎn)，比如rps12在分為三個區(qū)域且長度不一。因此诀浪，在合并不同提取結(jié)果中的同位點(diǎn)序列的時候還需要考慮拷貝問題。

   因此延都，作者也花了很多精力寫了一個R腳本（本文稱為多態(tài)性批量計(jì)算腳本）雷猪，可以針對曲小健腳本生成的多個提取結(jié)果，進(jìn)行同基因合并晰房、muscle多序列比對求摇、變異參數(shù)計(jì)算（包括序列長度、比對gap數(shù)量殊者、變異位點(diǎn)和變異比例与境、核苷酸多態(tài)性等）、基因名稱規(guī)范（針對perl腳本的命名風(fēng)格）的自動化實(shí)現(xiàn)猖吴，并生成一個csv匯總表格摔刁。本腳本還會同時輸出每個位點(diǎn)的比對好的fasta文件，放置在同路徑的fasta文件夾內(nèi)海蔽，可用于后續(xù)可能的合并分析共屈，比如針對類群關(guān)系較遠(yuǎn)的類群，僅用葉綠體基因組的cds區(qū)域或者genes來建樹党窜。腳本中默認(rèn)不計(jì)算比對長度小于200bp的位點(diǎn)拗引，也可以根據(jù)自己的需要對腳本進(jìn)行修改。

   具體使用上比較簡單幌衣，將多個樣品的perl腳本提取結(jié)果放到同一個文件夾中矾削，然后在R中或者cmd中使用腳本即可，會實(shí)時顯示運(yùn)算結(jié)果豁护。

# Date： 2022-08-05
# Author： Zhang Zhen
# Email： ficuszz@163.com
# Website: zhangzhen.plus
# Description： 讀取當(dāng)前文件夾的fasta文件哼凯，根據(jù)文件名提取出同源序列，批量比對楚里，計(jì)算變異位點(diǎn)比例挡逼，然后寫入scv文件
#               fasta文件準(zhǔn)備方式：使用曲小健博士的Perl腳本從葉綠體基因組gb文件中提取gene、CDS或者間隔區(qū)序列
#               gb文件準(zhǔn)備：下載葉綠體基因組全長序列腻豌，利用曲小健的PGA程序重新批量注釋基因組家坎，可保證所有注釋結(jié)果風(fēng)格類似，
#                          各基因和片段名亦相同

# 載入所需庫
library(stringr)
library(ape)
library(pegas)

# 第一部分吝梅，批量讀取fasta文件內(nèi)容
# 批量讀取當(dāng)前工作目錄中的fasta文件虱疏，然后再根據(jù)fasta文件名整理出物種名列表，最后將文件中的fasta文件逐個復(fù)制給物種名列表
# 此處注意fasta文件的命名苏携，必須統(tǒng)一為“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])
}
# 在當(dāng)前目錄下新建fasta文件夾做瞪，用于存儲后續(xù)的單基因比對文件
if (dir.exists("fasta")){
  unlink("fasta", recursive = TRUE)
  dir.create("fasta")
} else dir.create("fasta")





# 第二部分，指定一些序列名稱和序列合并的函數(shù)
# 處理fasta文件中的序列名，刪除perl腳本提取gene和intergenic時附帶的物種名装蓬，
# 使得來自不同文件的同源基因名字相同著拭，便于后續(xù)比較提取
seqname_deal <- function(x){
  str_sub(x,1,str_locate("_", string = x)[1]-1)
}

# 相反的處理方式，去除序列Labels包含的序列名牍帚，僅含物種名儡遮，后續(xù)批量構(gòu)建單基因樹的時候可保證tip labels一致，便于比較
seqname_deal1 <- function(x){
  str_sub(x,str_locate("_", string = x)[1]+1,nchar(x))
}

# 計(jì)算突變位點(diǎn)數(shù)量暗赶，變異位點(diǎn)和indel位點(diǎn)的綜合
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區(qū)域的重復(fù)片段，且包括rps12反式剪切體

# 處理rps12的特殊情況蹂随，rps12在葉綠體上有三個區(qū)域,其中一個較短十嘿，約140bp，剩余兩個對稱位于IR區(qū)岳锁，此處刪除掉較短的那個
# 避免對后續(xù)的多樣性分析結(jié)果造成影響绩衷，此處篩選的標(biāo)本是小于400bp，兼顧適當(dāng)?shù)拈L度變異
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] 
}    

# 處理重復(fù)序列唇聘，對所有的序列名稱內(nèi)的字符進(jìn)行重排，重排后一致的為重復(fù)序列
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))

# 根據(jù)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])
    }
  }
}
# 執(zhí)行刪除
if(length(delete.list) != 0) {
  name_list[[1]] <- name_list[[1]][unique(delete.list)]
}





# 第四部分柱搜，處理一些個人指定需要刪除的序列
# 定義需要刪除的結(jié)果迟郎，包括IR的反向重復(fù)基因，以及對含有內(nèi)含子的基因的處理,如：
# trnK,rps16,trnG,ycf3,trnL,trnV,clpP,rpl16,trnI,trnA編碼區(qū)過短聪蘸，長度主要由內(nèi)含子組成宪肖，在此不計(jì)入基因中
# ycf68為假基因，位于trnL的內(nèi)含子中健爬，在此僅計(jì)入trnL的內(nèi)含子控乾，不計(jì)算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]
}





# 第五部分，程序主體娜遵，計(jì)算各類參數(shù)如比對長度蜕衡，核苷酸多態(tài)性，gap數(shù)量
# 初始化
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()

# 執(zhí)行程序主體设拟，逐個計(jì)算同源序列的多樣性參數(shù)
for(i in 1:length(name_list[[1]])){
  temp_dna <- name_list[[1]][i]
  # 在第2至最后一個fasta文件中慨仿，找第1個fasta文件第i個基因片算的同源序列，并寫入變量中纳胧，待下一步計(jì)算
  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")
  # 逐個計(jì)算基因長度镰吆、變異位點(diǎn)數(shù)量、gap數(shù)量跑慕、總數(shù)量万皿、比例和核苷酸多態(tài)性
  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),
              "突變位點(diǎn)數(shù)量(含indel)為:",substi_number(align_temp_dna),
              "突變位點(diǎn)數(shù)量比例為:",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))






# 第六部分摧找，規(guī)則化輸出結(jié)果，包括列名和基因名
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")
# 替換不規(guī)則序列名
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]
}
# 寫入結(jié)果到文件
write.csv(last_result,"last_result.csv", quote = F,  row.names = F)
print("請?jiān)诋?dāng)前工作目錄下的last_result.csv文件中查看結(jié)果")

需聲明牢硅，此腳本仍非常稚嫩蹬耘，如有使用者，需要自行負(fù)責(zé)可靠性减余，可能也需要自行修改部分代碼综苔。

3. 完成做圖和可視化

   輸入結(jié)果稍作整理（格式如下），就可以使用以下腳本來完成佳励。

   
   

腳本如下：

library(ggplot2)
library(cowplot)
library(showtext)
install.packages("showtext")
font_add('Arial','c:/Windows/Fonts/arial.ttf') # 解決字體加載錯誤問題
showtext_auto() # 自動調(diào)用showtext休里，否則無法在ggsave()中使用蛆挫，因?yàn)間gsave會自動打開和關(guān)閉圖形設(shè)備


gene <- read.table("gene.txt",header = T,sep = "\t")
# 按照出現(xiàn)順選指定橫座標(biāo)順序赃承，預(yù)防下一步作圖是ggplot自行排列順序 http://www.reibang.com/p/437906c6b063
gene$lociname <- factor(gene$lociname,levels = unique(gene$lociname))

# 全部數(shù)據(jù)作圖
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位位點(diǎn)標(biāo)記獨(dú)特顏色
  geom_hline(yintercept = mean(gene$nucleotide.diversity.Pi)) + # 添加橫線悴侵，橫線位置為所有數(shù)值的平均值
  annotate("text",x = "rrn16", y = 0.0019,label = "0.001676",size = 4, hjust = 0.65) +  # 添加一個單獨(dú)的數(shù)字到圖中瞧剖，概述中為橫線的具體數(shù)值
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), text = element_text(size =14)) + #定義橫座標(biāo)字號和對其位置
  labs(x=NULL,y=NULL) +  # 不使用橫縱坐標(biāo)說明位子
  geom_text(label = gene$alignment.length,vjust = 0.4, hjust = -0.3, 
            position = position_dodge(1.5),size = 4, angle = 90) + # 條形圖上添加比對長度數(shù)值
  scale_y_continuous(expand = expansion(mult = c(0.02,0.12))) # 調(diào)整灰色畫布的上下邊界，利于一些標(biāo)注可免，使其不出框

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")


# 導(dǎo)出PDF文件
ggsave("test.pdf",width = 18,height = 10)

出圖效果如下抓于，包括基因和基因間隔區(qū)兩種分類、比對長度浇借、平均核苷酸多態(tài)性捉撮，以及多態(tài)性top5的藍(lán)色突出顯示。

2.7 選擇壓力分析

選擇壓力分析是作者之前基本沒有涉及的領(lǐng)域妇垢，聽說過的也就是CodeML巾遭、easycodeml以及datamonkey等，這些軟件的理論較為復(fù)雜闯估，而且也似乎不能批量的完成多位點(diǎn)多樣品的選擇壓力分析灼舍，在葉綠體基因組比較分析中較難應(yīng)用。而多個位點(diǎn)多樣品的選擇壓力分析本身也較為復(fù)雜涨薪，需要以下步驟完成：

1. 提取出注釋好cds序列

2. 合并多個樣品的同一cds位點(diǎn)的序列并比對

3. 針對任一個cds位點(diǎn)骑素，計(jì)算參考序列和所有樣品序列之間的Ka（非同義替換）和Ks（同義替換）兩個參數(shù)，通過計(jì)算Ka/Ks來評估選擇壓力

4. 批量完成所有cds位點(diǎn)的Ka/Ks的計(jì)算

5. 可視化

對于第一步和第二步刚夺，和多態(tài)性位點(diǎn)分析中類似献丑，同樣先使用曲小健的perl腳本提取所有注釋gb文件的cds序列，然后使用作者的多態(tài)性批量計(jì)算腳本進(jìn)行逐個位點(diǎn)的cds序列合并并比對侠姑。

對于第三步和第四步的任務(wù)阳距，需要考慮批量計(jì)算的可能行，作者選擇了TBtools的simple Ka/Ks Calculator功能完成Ka/Ks的計(jì)算结借。該功能需要兩個輸入文件筐摘，一個是比對好的fasta文件，一個是序列名對的二列的表格，用來指定計(jì)算對象咖熟，因?yàn)槊看蜬a/Ks的計(jì)算都是針對兩條序列的圃酵。注意，因?yàn)榇瞬叫枰獏⒖夹蛄锈晒埽栽诘谝徊胶偷诙侥抢锕停托枰獙⒖蓟蚪M一同運(yùn)算在內(nèi)。那么确沸，批量計(jì)算的步驟就很明顯了：

1. 合并所有cds位點(diǎn)的fasta文件捌锭，也就是將作者腳本生成多個比對好的fasta格式的cds序列文件合并到一起。

2. 生成一套Ka/Ks序列名稱對罗捎，用以指導(dǎo)TBtools的計(jì)算观谦。這套名稱對，應(yīng)該為每個cds位點(diǎn)桨菜，都指定好參考序列和所有的樣品序列之間名稱對豁状。這個名稱對的生成作者也利用R腳本生成，輸入數(shù)據(jù)即為作者多態(tài)性批量計(jì)算腳本輸出的fasta文件夾倒得。該腳本也會同時合并所有的cds序列泻红，可直接作為TBtools Ka/Ks計(jì)算的輸入文件。注意霞掺，此腳本要求參考序列位于所有樣品之后谊路。

   library(ape)
   # 批量獲取文件夾內(nèi)文件名
   name_list <- list.files()
   # 批量讀取文件夾內(nèi)的全部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)){
     # 給采集號為主的序列名加上位點(diǎn)名稱后綴
     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)}
     # 按位點(diǎn)逐個生成序列名稱對
     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為含參考序列的樣品總數(shù)，需自行調(diào)整
   }
   # 轉(zhuǎn)置列表形式的序列名稱對至矩陣形式
   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)
   

3. 將腳本生成序列對名稱文件和合并好的fasta文件使用TBtools來計(jì)算Ka/Ks即可菩彬。

第五步可視化缠劝，而可視化涉及到一個問題。一般來說挤巡，Ka/Ks大于1剩彬，被認(rèn)為是正向選擇；小于1矿卑，反向選擇或者純化選擇喉恋；等于1，中性選擇母廷。但是轻黑，在實(shí)際運(yùn)算中，還包括一些異常值琴昆，這些異常值在近緣類群中也不算少見氓鄙，但卻在已發(fā)表文獻(xiàn)中討論較少。比如业舍，Ka和Ks都是0的時候抖拦，軟件會給出NaN值升酣，這種一般也可以看作為中性選擇，畢竟Ka和Ks都沒有發(fā)生态罪；還有就是Ka大于0而Ks等于0的時候噩茄，會出現(xiàn)無窮大的結(jié)果，盡管也可以認(rèn)為是大于1的正向選擇复颈，但是在可視化中卻難以處理斗遏。

此外匕荸，在葉綠體基因組比較分析中，常常由于多個樣品工窍，每個位點(diǎn)會有多個Ka/Ks值箱玷。而文獻(xiàn)中通常僅展示一個值碌冶，可能為其平均值溺森，但求其平均也會降低信息含量癞己。又考慮到上段提到的異常值問題，作者給出一個較少見的箱形圖可視化方案舒帮，同時展示多個Ka/Ks值的分布特征和異常值会喝，具體如下陡叠。上圖為無窮大異常值的數(shù)量分布玩郊，下圖多個Ka/Ks值的箱形圖，其中NaN值被視為1枉阵。

該圖的生成需要準(zhǔn)備三個數(shù)據(jù)文件和一個做圖腳本译红。

數(shù)據(jù)文件1，剔除掉所有inf值的Ka/Ks結(jié)果表兴溜，NaN已被替換為1侦厚。該文件用來完成箱形圖。

數(shù)據(jù)文件2拙徽，記錄每個cds位點(diǎn)inf值出現(xiàn)次數(shù)的表格刨沦。該文件用來完成記錄inf值出現(xiàn)次數(shù)的折線圖。

數(shù)據(jù)文件3膘怕，是數(shù)據(jù)文件3的簡化想诅，僅記錄inf值出現(xiàn)次數(shù)不為0的cds位點(diǎn)的結(jié)果。該文件用來標(biāo)注折線圖具體代表的出現(xiàn)次數(shù)岛心，0次不標(biāo)注来破。如果你對R語言做圖非常熟悉，這個表格應(yīng)為多余忘古，可通過命令來完成0次不標(biāo)注的效果徘禁。但作者R水平有限，采取另準(zhǔn)備一個文件的笨辦法髓堪。

而所需做圖腳本如下：

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)) + # 調(diào)整橫座標(biāo)
  # 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) # 添加額外的備注（正無窮位點(diǎn)的值）
kaks

ggsave("kaks.pdf",width = 14,height = 5)

3. 后記和聲明

本文記錄所記錄的分析工具和可視化方案送朱，均為作者在發(fā)表文章中所用樣式娘荡，僅為個人偏好，讀者可酌情參考驶沼。

此外它改，盡管作者在本文中列舉了諸多個人計(jì)算和繪圖所用R腳本，但是其質(zhì)量商乎、效率央拖、可靠性和通用性都不敢有所保證，僅為作者數(shù)據(jù)分析時所用工具鹉戚。如需參考還需根據(jù)個人的工作路徑鲜戒、樣品數(shù)量或者測序文件性質(zhì)針對性的修改。所以可能仍需要讀者具有一定的R語言基礎(chǔ)抹凳，這些半成品腳本才有可參考的價值遏餐。因此，標(biāo)題亦有夸張之嫌赢底，僅為吸引讀者失都，還請見諒！