1.安裝samtools并將sam文件轉(zhuǎn)為bam文件

samtools是一款非常好用的工具，本文我們使用的功能較少。詳細(xì)的使用方法請見：samtools用法詳解_samtools 使用_sunchengquan的博客-CSDN博客

conda create -n samtools
#創(chuàng)建安裝環(huán)境
conda activate samtools
#激活環(huán)境
conda install samtools
#通過conda安裝samtools
cd /to/your/work/place
#轉(zhuǎn)到存放有sam文件的文件夾，sam文件為getorganelle組裝葉綠體基因組時(shí)產(chǎn)生，位于seed文件夾下
samtools view -@8 -b test.sam > test.bam
#使用samtools view將sam文件轉(zhuǎn)為bam文件,-@表示使用的線程數(shù)

2.排序歸檔和提取

# 構(gòu)建索引
samtools sort -o test_sort.bam test.bam
#對bam文件進(jìn)行排序
samtools index test_sort.bam test_sort.bam.bai
#對test_sort.bam構(gòu)建索引
samtools depth test_sort.bam > test_depth.txt
#計(jì)算每一個(gè)位點(diǎn)或者區(qū)域的測序深度并生成文件儲存
cut -f 2-3 test_depth.txt > depth.txt
#對生成的txt文件進(jìn)行數(shù)據(jù)裁剪，直接使用Linux的cut命令，-f表示裁剪整列，2-3表示裁剪后保存2-3列

3.合并步長，精簡信息

import os
# 腳本只需要修改file_path 和 output_file

# output_file是2000步長的合并文件
#file_path是上一步生成的深度文件depth.txt
file_path = open(r'D:/python/depth.txt', 'r')
output_file = r'D:/python/depth_2000.txt'

cont = file_path.readlines()

with open(output_file, 'w') as ff: # output file
    for i in range(len(cont)):
        if i % 2000 == 0:
            start = i
            end = min(i + 2000, len(cont))  # 確保不超過列表的長度
            middle = start + (end - start) // 2  # 計(jì)算實(shí)際區(qū)間中點(diǎn)
            all_sum = 0
            for j in range(start, end):
                all_sum += int(cont[j].split('\t')[2])
            average_depth = round(all_sum / (end - start), 0)
            ff.write(f"{middle}\t{average_depth}\n")

4.畫測序深度圖

先將上一步得到的count文件放入R的工作目錄

#查看工作目錄
getwd()
#設(shè)置工作目錄
setwd("/path/to/your/work/place")
#如果沒有安裝ggplot2，可以安裝install.packages(c("package1","package2",....)),此命令可同時(shí)安裝多個(gè)包
install.packages("ggplot2")
#加載包
library(ggplot2)

#加載數(shù)據(jù)
data <- read.delim("D:/python/depth_2000.txt", header=FALSE, sep = '\t')
data2 <- read.delim("D:/python/depth.txt", header=FALSE, sep = '\t')

##按照2000的步長間隔進(jìn)行繪圖
ggplot(data, aes(V1, V2)) + 
  geom_bar(stat = 'identity', width = 800, fill = "#BBE7FC") + 
  ylim(0,1500) + theme_classic() + xlab("Sequence length") + 
  ylab("Mean base depth")+
  geom_rect(aes(xmin=1000, xmax=86319, ymin=0, ymax=0.5), 
            fill="#BBE7FC", colour="#FDDBB2", linewidth=1.5) + #xmax=86319是LSC的終點(diǎn)
  geom_rect(aes(xmin=86320, xmax=111890, ymin=0, ymax=1), 
            fill="#B45218", colour="#B66300", linewidth=1.5) + #xmin=86320是IRB的起點(diǎn), xmax=111890是IRB的終點(diǎn)
  geom_rect(aes(xmin=111891, xmax=130299, ymin=0, ymax=0.5), 
            fill="#F3D68A", colour="#F4D68A", linewidth=1.5) + #xmin=111891是SSC的起點(diǎn), xmax=130299是SSC的終點(diǎn)
  geom_rect(aes(xmin=130300, xmax=155870, ymin=0, ymax=0.5), 
            fill="#C69D3C", colour="#C69B3C", linewidth=1.5) #xmin=130300是IRA的起點(diǎn), xmax=155870是IRA的終點(diǎn)


##按照所有位點(diǎn)繪圖
ggplot(data2, aes(V1, V2)) + 
  geom_bar(stat = 'identity', width = 800, fill = "#BBE7FC") + 
  ylim(0,1500) + theme_classic() + xlab("Sequence length") + 
  ylab("Mean base depth")+
  geom_rect(aes(xmin=1, xmax=86319, ymin=0, ymax=0.5), 
            fill="#BBE7FC", colour="#FDDBB2", linewidth=1.5) + #xmax=86319是LSC的終點(diǎn)
  geom_rect(aes(xmin=86320, xmax=111890, ymin=0, ymax=1), 
            fill="#B45218", colour="#B66300", linewidth=1.5) + #xmin=86320是IRB的起點(diǎn), xmax=111890是IRB的終點(diǎn)
  geom_rect(aes(xmin=111891, xmax=130299, ymin=0, ymax=0.5), 
            fill="#F3D68A", colour="#F4D68A", linewidth=1.5) + #xmin=111891是SSC的起點(diǎn), xmax=130299是SSC的終點(diǎn)
  geom_rect(aes(xmin=130300, xmax=155870, ymin=0, ymax=0.5), 
            fill="#C69D3C", colour="#C69B3C", linewidth=1.5) #xmin=130300是IRA的起點(diǎn), xmax=155870是IRA的終點(diǎn)

Rstudio可以選擇輸出文件的類型，有png和pdf兩種。

完
本文參考鏈接為https://blog.csdn.net/salty_fish_xu/article/details/136592818?spm=1001.2014.3001.5502